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Introduction. The problem of irrotational gravity waves in water is, from the 
mathematical point of view, a problem in potential theory which involves a nonlinear 
boundary condition at the free surface. In addition, the shape of the free surface is 
itself not given a priori but is to be determined as part of the solution along with other 
quantities, such as the distribution of the velocity and the pressure. Very few success- 
ful attacks on the problem in this formulation have been made; among these are the 
proofs of the existence of steady periodic waves of special type by Levi-Civita [5 ],** 
and Struik [11], and a partial treatment of the problem of the so-called solitary 
wave by A. Weinstein [13]. 

Most of the literature on the subject of surface gravity waves is concerned with 
theories which result from the general nonlinear theory when simplifying assump- 
tions of one kind or another are made. The present paper is concerned with two such 
theories: 

1) Perhaps the best known and most extensively studied theory is that which re- 
sults from the general theory when it is assumed that the amplitude of the waves at 
the surface and the velocity of the particles there are small enough that the free sur- 
face condition can be simplified by dropping the nonlinear terms; in addition, this 
condition may be prescribed at the original undisturbed surface of the water. The re- 
sult is a problem in potential theory with a linear boundary condition of the mixed 
type. The greater part of our work here makes use of this theory. We shall refer to it 
as the exact linear theory, or simply as the exact theory. 

2) The second theory furnishes an approximation to the exact linear theory which 
is based on the assumption that the depth of the water is small. The resulting theory 
yields a differential equation for the surface elevation of the water which turns out to 
be the wave equation. This approximate theory, which is often called the shallow 
water theory is used, for example, in discussing the tides in the ocean. In Sec. 7 we 
give a brief derivation of the shallow water theory which brings out the role of the 
depth of the water as the determining factor in the accuracy of the approximation. 
The usual derivation of the theory based on assuming that the pressure in the water 
is given by the same law as in hydrostatics does not bring this point out clearly. The 
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approach to the theory through the hydrostatic pressure relation also does not lend 
itself easily to generalization to other cases, such as the derivation of the shallow 
water theory when floating bodies are present. 

Our principal object in this paper is to solve the problem of determining progress- 
ing waves over a uniformly sloping bottom by making use of the exact linear theory, 
to discuss the solutions numerically, and to compare them with the approximate solu- 
tions furnished by the shallow water theory. 

Solutions for waves on sloping beaches in terms of the exact theory have been ob- 
tained by Hanson [3], Bondi [2], Miche [9] and H. Lewy [6]. The first author ob- 
tained one type of standing wave solution. The second and third authors obtained 
two types of standing wave solutions for the case of motion in two dimensions from 
which progressing wave solutions can be constructed. The first three writers, as well 
as the writer of the present paper, are concerned only with cases in which the bottom 
slopes at the special angles 7/2n, with m an integer. The method employed in the 
present paper is different from those of the first three authors; in particular, the pro- 
gressing wave solutions are obtained here in a closed form which lends itself well to 
detailed discussion. Also, the method employed in the present paper yields three- 
dimensional progressing wave solutions, that is progressing waves which approach an 
arbitrary plane wave at infinity (cf. Sec. 9). 

The investigations which led to the present paper were begun in collaboration 
with H. Lewy, who then later extended the method to two-dimensional motions (cf. 
[6]) for a bottom sloping at the angles px/2n in which p is any odd integer and n is 
any integer such that 2n>p. (The cases p #1 are very much more complicated than 
those for p=1, by the way.) 

The basic idea of the method devised by H. Lewy and used by him and the author 
is to obtain a differential equation for the desired velocity potential, through use of 
the boundary conditions, which is not the potential equation and which, as it turns 
out, permits an explicit integration. In the case of two-dimensional motions the prob- 
lems can be treated, of course, by making use of analytic functions of a complex 
variable; in this case Lewy’s differential equation becomes an ordinary non-homo- 
geneous differential equation with constant coefficients for the complex potential, and 
this equation can be integrated to yield the desired solutions. In the three dimensional 
cases the solutions can also be obtained, as explained in Sec. 9, but the results are 
more complicated and more difficult to handle numerically. It may be of interest to 
observe that the method developed here for the mixed boundary value problem in 
wedge-shaped regions of angle 7/2 is not confined in its usefulness to solutions of 
the potential equation—it could also be extended to other linear partial differential 
equations. 

In Sec. 1 the exact linear theory (which apparently goes back to Poisson) is formu- 
lated briefly. In Sec. 2 the theory is applied to yield the well known solutions for 
steady progressing waves in water of infinite depth. In addition, we show in Sec. 2 
that these solutions are uniquely determined if the amplitude and velocity of the 
waves is bounded at ~.! The method used to obtain this result is essentially the same 
as that employed later to obtain progressing wave solutions over a sloping bottom. 

1 The author was led to these rather general conditions guaranteeing the uniqueness of the solutions 


as the result of a conversation with A. Weinstein, who had previously obtained the same result by quite 
different methods for water of finite and constant depth (cf. [12]). 
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In Sec. 3 convenient dimensionless independent variables are introduced, and 
these variables are used in the remainder of the paper except in Sections 7 and 8. 

In Sec. 4 the case of progressing waves coming from infinity in an ocean of infinite 
depth but bounded on one side by a vertical cliff is treated in considerable detail. 
Only the two dimensional case is considered—that is, the wave crests and all other 
curves of constant phase are assumed to be horizontal straight lines parallel to the 
cliff. Since the problem is then a potential problem in two dimensions, it is convenient 
to solve it in terms of analytic functions of a complex variable. In Sec. 5 the method 
used for the case of a vertical cliff is generalized to yield solutions in water over a 
plane bottom sloping at any angle 7/2n, with m an integer. Again only the two di- 
mensional case is treated. The essential step in the generalization requires the deriva- 
tion of Lewy’s differential equation, which turns out to be a differential equation of 
order m for an angle r/2n. The solutions of the differential equation which satisfy the 
boundary and regularity conditions are then given in this section. The solutions ob- 
tained are shown to be uniquely determined if they have at most a logarithmic 
singularity at the origin (that is, at the shore line) and satisfy certain boundedness 
conditions at ©.* The method of determining the arbitrary constants in such a way 
that the boundary conditions are satisfied is discussed in Appendix I. In Appendix II, 
the behavior of the solutions at © is investigated; we find that the solutions behave 
as one would expect, i.e. that they tend to the simple standing wave solutions for 
water of infinite depth. It is then readily seen that standing waves of arbitrary phase 
and amplitude at © can be constructed from our solutions. In this way the existence 
of a unique set of standing wave solutions is established. For the angles 7/2n, then, 
the problem of progressing waves has been solved completely for the two dimensional 
case. 

In Sec. 6 the theory of Sections 4 and 5 is applied to the cases n=1, n=2, and 
n=15, i.e. to the cases of bottom slopes of 90°, 45°, and 6° respectively.* The standing 
wave solutions are given numerically in the form of graphs for a distance of a few 
wave lengths from shore. The numerical evaluation requires the calculation of the 
values of complex integrals of the form 


+ 
E(z) = ef te—* dt 


taken over appropritate paths in the ¢-plane. A table of values of E(z) for the range of 
values of interest to us was computed and is included here as an Appendix. This table 
was based on a previous table calculated by the Mathematical Tables Project [8]. 
It might be noted that the calculations for the case of the 6° slope were very laborious. 

In all cases (except that of water of infinite depth everywhere) there are two types 
of standing waves: the one type has a finite amplitude all the way to shore, the other 
type has amplitudes which become logarithmically infinite as the shore is approached. 
At infinity the amplitude and wave length may be arbitrarily prescribed for both 
types, but the wave length at » and the frequency are connected by the same reta- 
tion as in water which is infinite in depth everywhere. These two types of standing 

2 This result (which is not obtained in the papers by Miche, Bondi, and Lewy cited above) means 
that no non-trivial solution exists which dies out at ©. 

* The paper of Miche [9], which gives the solutions in the general case, contains graphs for the 45° 
and 30° cases as well as approximate solutions for small angles of slope. 
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waves are always out of phase at ~, but the ones remaining finite (for any fixed bot- 
tom slope) at the shore all have the same phase at infinity. (The solutions obtained 
by Hanson [3] are those which remain finite at the shore.) Thus all progressing waves 
furnished by the exact theory necessarily have large amplitudes near shore. 

In all three cases treated numerically, i.e. those with 90°, 45°, and 6° slopes, the 
most striking general result is the following: The wave lengths and amplitudes change 
very little from the values at © until points about a wave length from shore (wave 
length at © is meant) are reached. Closer in shore the amplitude of any progressing 
wave becomes large. It is curious, however, that the amplitude (that is, the relative 
maximum or minimum) of a progressing wave at a point a wave length or two from 
shore can actually be about 10 per cent Jess than the value at ©. This was found in 
all three cases, including the 6° case. 

In Sec. 7 the shallow water theory is derived. In this theory, the wave amplitudes 
for a uniform bottom slope die out at © like x—'/4 (where x denotes distance from 
shore) while the wave lengths (or, rather, the distances between successive nodes) 
increase without limit. Thus if the wave lengths and amplitudes were to be taken with 
the same value at a given point some distance off shore in the two theories (i.e. the 
exact theory and the shallow water theory) the amplitude near shore as given by the 
shallow water theory would be much too high while at great distances from shore it 
would be too low. 

In Sec. 8 the results of the exact theory are compared with those of the shallow 
water theory for the case of water of finite and constant depth in order to bring out 
the known fact that the accuracy of the shallow water theory for a simple progressing 
wave depends upon having the ratio of depth to wave length sufficiently small. 

Progressing waves moving toward shore as given by both theories for the case 
of the 6° beach are compared graphically in Sec. 8 assuming the same frequency in 
both cases and also the same phase and amplitude at a point two or three wave 
lengths from shore.The results from there on toward shore do not differ greatly ex- 
cept at points very close to shore. This was to be expected, since the depth of the water 
at the point where the phase and amplitude are the same is about one-eighth of the 
wave length and hence one might expect the shallow water theory to be quite accu- 
rate. However, if the amplitudes had been made equal at a point 15 or 20 wave lengths 
from shore, the amplitudes given by the shallow water theory three wave lengths 
from shore would be 50 per cent to 60 per cent higher than those furnished by the 
exact theory. In other words, the amplitude variation with decrease in depth cannot 
be correctly estimated over distances of many wave lengths from a given point using 
the inverse fourth root law unless the wave length at the point is eight or ten times the 
depth of the water. On the other hand, as the results of the exact theory indicate, if 
the depth remains an appreciable fraction of the wave length the amplitude changes 
very little with changes in depth. We draw this conclusion from the fact that the wave 
amplitudes as given by the exact theory approach their values at © very quickly when 
the depth approaches say half the wave length. 

Another feature of our numerical results which is of interest concerns the variation 
in the phase or propagation velocity c of a progressing wave when the depth A varies. 
In Sec. 8 graphs are given which show the results of calculations of the phase velocity 
for the 6° case, using both the exact theory and the shallow water theory. In the 
case of the shallow water theory, it is found that the phase velocity c is practically 
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identical with that given by the formula c = (gh)'/*. In the case of the exact theory, the 
phase velocity is given very accurately by the formula c= [(Ag/2m) tanh (2xh/X) }1/2, 
that is, in both cases ¢ is given by the formula which is exact (for that theory) only 
when the depth / is constant. (In the second formula \ is taken as the wave length 
at «.) In other words, the phase velocity at any point is given very accurately by 
the formula which is exact only for water of uniform depth equal to that at the point 
in question and for a steady progressing wave traveling in one direction only. At the 
same time, the results for the phase velocity as furnished by the two theories agree 
very well with each other for a distance of about three wave lengths from shore but 
from then on out the shallow water theory furnishes a value which is too high by an 
amount which increases without limit with increasing distance from shore. These re- 
marks, it should be recalled, result from our calculations for a slope of 6°. With de- 
creasing slope, it seems certain that the shallow water theory would be accurate up 
to distances comprising more and more wave lengths from shore. 

All progressing wave solutions discussed above were obtained on the assumption 
that the wave comes from » toward shore with no component which goes outward 
at ©. Once the frequency and amplitude at © are prescribed, the additional condi- 
tion that the wave at © is a progressing wave moving toward shore leads to a unique 
solution, which, as we have already mentioned, has a logarithmic singularity at the 
shore line. The solution is also uniquely determined if the singularity at the shore line 
is prescribed—the behavior at © is then determined. Our theory thus furnishes us 
with two types of standing wave solutions from which solutions behaving like arbi- 
trary simple harmonic progressing waves at © can be constructed, but it furnishes no 
criterion by which one can decide what type of wave would actually occur in practice. 
Our assumption, in the numerical cases treated, that the waves move from © toward 
shore with no reflection from the shore back to © was the result of information on 
the phase velocities as measured on beaches with small slopes; these measurements 
agree rather well with the theoretical results discussed in the preceding paragraph 
for a progressing wave moving toward shore. The physical mechanism which prevents 
the reflection of waves from the shore can be understood as the result of the partial 
loss of energy from turbulence and the conversion of the remainder into an undertow 
through the occurrence of breakers. If, however, the slope of the beach is large it may 
well be that a standing wave, denoting perfect reflection, could occur. 

In Sec. 9 we solve the problem of progressing waves in an ocean of infinite depth 
bounded on one side by a vertical cliff when the wave crests are not assumed to be 
parallel to the shore line (as in Sec. 4); that is, we solve a three-dimensional problem 
using the exact theory. Solutions are obtained which tend at to an arbitrary plane 
wave. In all of the solutions obtained in the preceding sections by means of the exact 
theory the discussion was greatly facilitated by the use of analytic functions of a 
complex variable. In the present three-dimensional case this approach is no longer 
possible. Nevertheless, the process of obtaining the solutions remains analogous to 
that using complex functions. The solution for the case of a vertical cliff only is ob- 
tained, but it is readily seen how three-dimensional progressing wave solutions for 
slopes of angles 7/2n could also be found. The solution for the case of the vertical cliff 
is also evaluated numerically in Sec. 9 for the case of a progressing wave with wave 
crests tending to a straight line at © which makes an angle of 30° with the shore line. 
One of the figures given in Sec. 9 shows the contours of the wave surface. In this case, 
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as well as in the previous two-dimensional cases, it turns out that there is a point 
near the cliff where the wave crests are lower than they are at ©, although the eleva- 
tion of the wave crests becomes infinite upon approaching the cliff. 

Finally, the author takes pleasure in acknowledging the help and advice he re- 
ceived from a number of his colleagues and co-workers. The actual solution of Lewy’s 
differential equation and the determination of the constants to satisfy the boundary 
conditions—-no small task in itself—was carried out by E. Bromberg and E. Isaacson. 
The extensive numerical computations were completed by E. Isaacson, B. Gross- 
mann, and J. Butler. 

1. Resumé of general theory of surface waves of small amplitude. In this section 
we state briefly the well-known mathematical formulation of the problem of surface 
waves of small amplitude in water. (See, for example, Lamb: Hydrodynamics, Chap 
IX; or Milne-Thompson: Theoretical Hydrodynamics, Chap. XIV.) The water is as- 
sumed to fill the region —h(x, 2) SyS0 when at rest. The non-negative quantity h 
is the (variable) depth of the water. The motion is assumed to be irrotational, so that 
a velocity potential ®(x, y, 2; t) exists, in which ® depends not only on x, y, z but also 
on the time ¢. Hence ® satisfies the Laplace equation 
eb bh Bh 

+-—— = @, (1.1) 


Vb = 
Ox? ay? Oz? 








A solution of this differential equation in the region —hSy3S0 is to be found which 
satisfies appropriate boundary conditions. The condition to be satisfied at the free 
surface y=0 is 
0° OP il 
—-+g—=0, (1.2) 
Ot? oy 
which results from the Bernoulli law and the assumption that nonlinear terms in the 
displacement and velocity of the free surface can be neglected. At the bottom 


y= —h(x, 2) we require of course that the derivative of ® normal to the bottom sur- 
face should vanish: 
O® 
—- = 0. (1.3) 
On y=—h 


We shall be interested solely in phenomena which are periodic in the time ¢. It 
is therefore convenient to replace ® in the above equation by e**‘g, in which g depends 
only on x, y, and not ont. 

The conditions on ¢g are the same as those on ® except that (1.2) becomes 

Ya (4.2) 
oy 

In addition to the differential equation (1.1) and the boundary conditions (1.2’) 
and (1.3), the functions g(x, y, z) should be required to satisfy certain conditions at 
infinity which lead to unique solutions of the type desired on physical grounds. These 
conditions will be formulated in the following sections for the special cases of interest 
to us. 

Once ¢ has been determined, the vertical displacement n(x, 2; t) of the free surface 
is determined (see the books of Lamb and Milne-Thompson cited above) by the for- 


mula 
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1 0a 
n»=— —(e*'y); . (1.4) 
g ot ly=0 

2. Plane traveling waves in water of infinite depth. We seek solutions of the 
boundary value problem partially formulated in equations (1.1), (1.2’) and (1.3) 
which have the form of plane traveling waves for the case of infinite depth of water. 
The potential function g may be assumed to depend only upon x and y and not on z: 
yg =(x, y), and it is to be determined in the entire half plane y SO. 

The functions 


yg = e™ cos (mx + a), with m and a arbitrary, (2.1) 


are familiar potential functions which obviously yield “standing waves” upon rein- 
troduction of the time factor.4 The amplitude of these waves decreases exponentially 
with the depth. The condition (1.2’) at the free surface is satisfied if the following 


relation holds: 
2r 
lead nial, i“ (2.2) 


in which \ is the wave length. This yields, of course, a relation between the frequency 
and wave length which characterizes the type of dispersion encountered with surface 
gravity waves in water of infinite depth. 

Since our problem is linear and homogeneous, we may take linear combinations 
of the standing wave solutions to obtain “traveling wave” solutions given by the 
velocity potential ®(x, y, z, t) as follows: 


& = e™ cos (mx + ot + a), (2.3) 


in which @ is an arbitrary phase shift. This would represent a wave traveling in the 
direction of the negative x-axis. Of course, the relation (2.2) between wave length 
and frequency must always be satisfied. The above theory is the well-known classical 
theory, which is due to Poisson. | 

A point which seems not to have been raised in the standard treatises is that con- 
cerning the uniqueness of the solutions given by (2.1). It is of interest to deal with 
this question here since the reasoning used is later on generalized in such a way as to 
yield the solutions for the problem of waves in water over a uniformly sloping bottom. 

We wish to show the following: If g is a regular potential function in the half- 
plane y <0 which satisfies the free surface condition o*g = gdg/dy on y=0 and if the 
function g(x, y) and its first derivatives g, and g, are uniformly bounded in (x, y) 
as x?+y2— ©, then g=Ae”™” cos (mx+a), that is, if the velocity and the vertical dis- 
placement of the water are bounded at «, then ¢(x, y) is either identically zero or it is 
of the type (2.1) with A #0. The mixed boundary condition at the free surface in con- 
junction with the relatively weak condition at © thus leads to this rather narrow class 
of solutions.® 

Our proof of this theorm requires the introduction of the analytic function 


‘It is a curious fact that this theory, which deals with what is perhaps the most familiar of all wave 
motions in nature, is governed by the potential equation rather than the wave equation. Nevertheless, 
we shall use the terms wave-length, amplitude, phase, etc. with a meaning which is obvious from the context. 

5 The author’s attention was called to this possibility by A. Weinstein, who had already proved the 
same theorem for the case in which the fluid has a finite and constant depth. See [12]. 
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f(z) =¢+w of the complex variable® z=x-+1y, the real part of which is the potential 
function ¢. It is convenient at this point to translate our conditions at © on g(x, y) 
into conditions on f(z) at «. In both cases, of course, these functions are defined in 
the half-plane y <0. On account of the Cauchy-Riemann equations it follows at once 
that | df/dz| is not greater than | p2| + | $y and thus is uniformly bounded at © since 
Py| . Since f(z) = f#f’(£)dé it is clear that we have 


| f(z)| < M| 2], 
i 


where M a positive constant, for all sufficiently large ||. 

In addition, it is readily seen that | df/dz| tends to zero when |z|—>« along any 
ray in ySO0 which is not parallel to the real axis. Since we wish to make use of a 
slightly more general result later on it is convenient to formulate here the following 
Lemma: If f(z) =¢+i is analytic and regular in the interior of a sector of the complex 
plane and ¢ is uniformly bounded at ~, the absolute value of any derivative of f(z) tends 
to zero when | z| — 2 along any ray which 1s not parallel to either of the boundary rays 
of the sector. The lemma is an almost immediate consequence of the assumption that ¢ 
is uniformly bounded at «. We consider g(x, y) to be expressed in terms of its values 
on the circumference of a circle with center at (x, y) and radius R through use of the 
Poisson integral formula. By differentiating both sides of this relation one easily ob- 
tains bounds for ¢, and ¢, at the point (x, y) of the form | ge <2M,/R, | oy | <2Mi/R 
with M, a constant which may be taken as the maximum of || on the circle of 
radius R. It therefore follows that | df/dz| s2+iy<4Mi/R in view of the Cauchy- 
Riemann equations. We now assume that M, is a fixed upper bound for ¢ for all suffi- 
ciently large z—that such a bound exists was assumed. As | z| —« along any ray not 
parallel to the sides of our sector, it is clear that | df; dz|—>0 since we may allow R 
to tend to infinity (i.e. our domain will accommodate circles of arbitrarily large radius 
with centers on such a ray) while M, is fixed. Since g, and gy are potential functions 
which are well known to be bounded at ~ in any closed sub-domain of the sector’ 
it follows that the second derivatives of g tend to zero and therefore also | d*f/dz?| 
tends to zero along the rays in question. That the higher derivatives of f(z) behave in 





this is the case for |y,| and 





the same manner is now obvious. 
The condition on ¢ at the free surface leads to the following condition on f(z) 
=g+wy on the real axis. We may write successively 


a a 
(< —— *)¢ = Me (« —— °*) (g + ip) 
oy oy 
#4 
= te («i —— o*) f() 
dz 


(Re means “the real part of” the expression which follows) the last step being a con- 
sequence of the fact that f(z) is analytic. Thus the free surface condition may be 


written 


o 
a 


d 
Re («i —— o*) J = 0 for z real. (2.4) 
6 In all but the final section of this paper we deal with two-dimensional problems only so that no con- 
fusion between the complex variable z and the third space variable should arise, particularly since the 
complex variable is not used in the final section. 
7 The bounds 2M,/R for | ps| and | ¢,| obtained above could be used to yield this result. 
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We next introduce the analytic function F(z) defined by 





r d 
F(z) = («i —- o*) fe) (2.5) 
dz 
Clearly F(z) is regular in the lower half plane. Since f(z) satisfies (2.4) it follows that 
































the real part of F(z) is zero on the real axis, and hence F(z) can be continued analyti- 
cally by reflection into the entire upper half-plane; thus the resulting function is regu- 
lar in the entire plane, and any bounds for | F(z) | in the lower half plane also hold in 
the upper half plane. From (2.5) and our conditions on f(z) at © it is clear that 

F(z)| <Me|2\| for all sufficiently large |z| in the lower half plane, and hence also in 
the entire plane. It follows that F(z) is linear, by application of Liouville’s theorem to 
the function | F(z) — F(0) |/z; that is, it can be written as F(z) =cz+c2. If we now in- 
troduce c+ for F(z) in (2.5) and integrate to obtain f(z) the result is 





| f(s) = Ae™ + Bz +C, (2.6) 


with m=o*/g. The constant A is an arbitrary complex constant. The constant B 
must, however, be zero since |df/dz| tends to zero along all rays not parallel to the 
real axis by the Lemma proved above,’ and the same is true of | de- im? /dz| since m 
is real and §te( —imsz)—+— © along such rays. The constant C must be pure imaginary 
because of the boundary condition (2.4), as one readily sees. Thus the only non-van- 
ishing potential functions g for which the velocity and surface elevation are bounded 
at « are of the form (2.1). 

3. Introduction of dimensionless quantities. In dealing with surface waves in the 
remainder of this paper it is convenient to work with dimensionless space variables 
x, and y; defined by x; =mx, y:=my, in which m is given by 

2r 
o? = gm = g-—» (3.1) 
r 
with o the circular frequency and \ as wave length. We also replace the time variable ¢ 
by a new variable 4; given by t;=ot. In these variables the surface condition (1.2’) 
is readily seen to take the form 








deg 
—-—g=0 for y=0, (3.2) 
Ovi 
where ¢ of course satisfies V?9 =0 in the new variables. However, since nearly all of 
our work from now on is carried out in the new variables, we shall drop the subscripts 
but retain the surface condition in the form (3.2). The original variables can always be 
reintroduced by replacing x and y in all of our results by mx and my and t by at. 
This means that the reciprocal of the quantity m defined by (3.1) is our unit of length. 
In the course of our discussion on waves over a sloping bottom it will be shown that 
the relation (cf. (2.2)) ¢2=gm=2rg/X between frequency and wave length for water 
of infinite depth holds asymptotically as the depth of the water becomes infinite, 
when d is the “wave length at «.” Thus our unit of length in these cases is propor- 


tional to the wave length at «. 





8 At this point we use the assumption that ¢ is bounded. 
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The standing wave solutions ® corresponding to (2.1) are, in the new variables 
(after dropping subscripts) : 
@ = e''ev cos (x + a). (3.3) 


4. Plane traveling waves in an ocean of infinite depth bounded on one side by a 
vertical cliff. As stated in the introduction, the main purpose of this paper is to study 
plane traveling waves in an ocean with a uniformly sloping bottom. In this section 
we deal in detail with the special case in which the “bottom” is vertical. Most of the 
essentials of the method to be employed in the more general cases are well illus- 
trated in this case, while the formal apparatus necessary is much simpler than that 
needed for the general case. In our treatment of the more general cases we shall then 
feel free to condense the presentation in many particulars. 

We assume that all quantities depend upon x and y only so that all curves of con- 
stant phase (the loci of the wave crests, for example) are parallel to the line of inter- 
section of the free surface and the cliff forming the vertical boundary of the water. 
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Thus we seek a potential function g(x, y) in the shaded area of Fig. 1 which satisfies 
the surface condition (see Sec. 3): 


00 
—-—g=0 when y=0, x>0 (4.1) 
Oy 
and the condition at the vertical wall 
00 
a = 0 when x=0, y <0. (4.2) 
x 


As we have already stated, our purpose is to obtain potential functions g which 
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satisfy (4.1) and (4.2) and which behave at like traveling waves moving toward 
shore. It seems reasonable to expect that a velocity potential ®(x, y, ) which behaves 
at © like the solution ee” cos (x +a) for water of infinite depth everywhere (cf. (3.3)) 
will exist in the present case. Or, as we could also put it, we may expect that two po- 
tential functions ¢;, g2 can be found which behave on the surface at © like sin x and 
cos x respectively, since a reintroduction of the time factor e* would yield two “stand- 
ing wave solutions” which could be combined linearly (since our problem is linear and 
homogeneous) to yield a solution behaving like a traveling wave at .° In what fol- 
lows we shall obtain two such potential functions which are 90° “out of phase” at ©. 
One such solution which is bounded and regular can be obtained immediately: The 
boundary condition (4.2) permits an analytic continuation of ¢ by feflection in the 
negative y axis into the entire half plane y £0, and we have already obtained solutions 
for this case in Sec. 2. Since only an even function of x is in question it follows that 
¢:1 = Ae” cos x is the only solution regular in the entire fourth quadrant which, together 
with its first derivatives, is uniformly bounded at © in this quadrant, because of the 
fact that the solutions obtained in Sec. 2 were shown to be unique under these 
circumstances. (It is clear that bounds for g at © in the fourth quadrant hold, after 
reflection, in the half plane y $0.) In other words, all non-singular solutions which are 
bounded at ~ have the same phase at ~, i.e. they behave like cos x there.'® To obtain 
solutions “out of phase” with e” cos x at © it is therefore essential to admit a singu- 
larity. On the other hand, it is rather natural on physical grounds to expect a singular- 
ity at the origin (i.e. at the water line on the vertical cliff) of the type of a source or 
sink if a progressing wave coming toward shore from © occurs.' This point has al- 
ready been discussed in the introduction. 

We are now in a position to complete the formulation of our boundary value prob- 
lem by prescribing conditions on ¢ at the origin and at ~. We require, in accordance 
with the remarks above, that g should possess a representation of the form 


o= ¢ log rT + re (4.30) 


valid near the origin, with r=x?+y? and g and g functions which together with their 
first two derivatives are bounded in a neighborhood of the origin. At © we require 
that ¢ and its partial derivatives of the first two orders be uniformly bounded, i.e. 


that a constant M exists such that 
lel + Dleo™| < M, (4.3,.) 


for all sufficiently large x?+y?, in which the sum is taken over all first and second 
partial derivatives of g.'* The conditions (4.1), (4.2), (4.39) and (4.3,,) constitute the 


® Of course the motion will not be a steady wave motion in general, but one which “approaches” 
a steady motion at ~, 

10 Upon reintroduction of the original space variables it is seen that this type of solution includes 
waves of all possible wave lengths. 

11 The paper of Hanson [3] mentioned in the introduction contains only the regular solution. It 
might be of interest to note that the starting point of the present investigation was the conjecture that 
solutions out of phase with those of Hanson at © could be obtained by admitting a singularity corre- 
sponding to a source or sink at the origin. 

12 These requirements are more stringent than would be necessary to ensure the existence and unique- 
ness of the type of solutions desired. However, we are not interested in this paper in formulating condi- 
tions at © and at the origin in the most general way possible, but only in formulating conditions which 
seem reasonable on physical grounds and which will lead to unique solutions of a type which interest us. 
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complete set of conditions on g. We shall obtain all non-vanishing solutions of this 
problem by constructing them explicitly. . 

Our method of solving this boundary value problem requires the introduction of 
‘the analytic function f(z) of the complex variable z=x+iy, the real part of which is 
the potential function ¢: 

f(z) =e + w. 
We must reformulate our conditions on ¢ in terms of f(z). The boundary conditions 
(4.1) and (4.2) can be written as 


0 fe) E 
( — i) = te (— ~ 1) + ww) 
Ov Oy 


d 
= ate (4 —- 1) A = ReL,(D)f = 0 on the positive real axis (4.4) 


dz 
and 

0 0 2s 

—(¢) = Re— (e+ wy) 

Ox Ox 

df(z) ee 
= Re - = ReL.(D)f = 0 on the negative imaginary axis. (4.5) 

az 


The last step in each case is justified by the fact that f(z) is analytic. The symbol Re 
means, of course, that the real part of what follows is to be taken, and the symbols 
Li and Lz refer to the linear operators defined by (4.4) and (4.5). (D=d /dz.) 

The conditions (4.3,,.) at © on ¢ lead to the following conditions on f(z) at «: 
1) Because of the fact that ¢ is uniformly bounded at o, the Lemma of Sec. 2 
shows that | df /dz| and |d*f/dz?| tend to zero along all straight lines in the quarter- 
plane which are not parallel to its boundaries. 2) Because of the Cauchy-Riemann 
equations | df/dz| and | d°f/dz?| are uniformly bounded in z as | z| 0, 

We shall make use of the condition (4.39) at the origin in the following form: The 
analytic function f(z) is such that | df /dz| <M,/|z and | d2f/dz*| <M; | z| 2 with M, 
and M; constants, in a neighborhood of the origin. This statement follows immedi- 
ately from the conditions (4.39) since r¢@z, rey, 7°22, 1*Gzy are as a result bounded near 
the origin and this leads to finite bounds for |zf’| and |s*f’’| through use of the 
Cauchy-Riemann equations. 

Our method of solution depends essentially upon the observation that the linear 
operators L; and Lz defined by (4.4) and (4.5) have the following property : 





MeL, L2(f) = ReLoL,(f) = (0) (4.6) 


on both boundaries of our domain, i.e. on both the positive x-axis and the negative 
y-axis. This property is an immediate consequence of the linearity and special form 
of ZL; and Lz and the boundary conditions (4.4) and (4.5). We are thus led to consider 
the analytic function F(z) defined by (see (4.4) and (4.5)) 


d d 
F(z) = LoLlif(z) = (i — 1) f(). (4.7) 


dz 
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We shall prove the following: If the analytic function f(z) having the properties we have 
postulated exists, then F(z) =Ai/z*, with A an arbitrary real constant; hence f(z) would 
satisfy an ordinary differential equation with constant coefficients. 

What are the properties of F(z), assuming that our f(z) exists? From (4.6) we ob- 
serve that the real part of F(z) vanishes on both boundaries of the quarter-plane (i.e. 
on the positive real axis and the negative imaginary axis). F(z) can therefore be con- 
tinued analytically by the reflection process into the entire plane; the result will ob- 
viously be a single-valued function whose real part vanishes on the entire real as 
well as the entire imaginary axis. (Here we make essential use of the fact that the 
original domain is a sector of angle 7/2.) At the origin F(z) has at most a pole of 
order two since | df /dz| <M,/|z| and |d*f/dz? <M-2/|z| 2 hold near z=0, and these 
bounds for the derivatives of f in the quarter-plane lead to the bound | F(z)| < M,/|:| . 
in a full neighborhood of the origin as one sees from (4.7) and the fact that F(z) is 
continued by reflection into a single-valued function in the entire plane. Hence F(z) 
has at most a pole of order two, and not an essential singularity, at z=0. In the same 
way the conditions at » on f(z) yield for F(z) the condition that | F(z)| is uniformly 
bounded at «. Also | F(z)| tends to zero when |z|—>« along any ray which is not 
parallel to the real or the imaginary axis, since | df/dz| and | d?f/dz?| have this prop- 
erty. The only analytic function F(z) with all of these properties is F(z) = Ai/z*, with 
A an arbitrary real constant (zero included): It is well known that a single-valued 
analytic function defined in the entire plane is determined by its singularities, which 
in this case consist of a pole at the origin. This fact, together with the additional con- 
ditions on F(z), leads easily to our result. 

We can now be certain that the solutions g(x, y) of our potential problem must be 
the real part of an analytic function f(z) which satisfies the ordinary differential equa- 











tion 


d r d 1 
- (: — 1) A) =A—, A real. (4.8) 
dz\ dz , 


oe 
. 


Our problem is therefore reduced to that of determining integration constants 
in the solution of (4.8) in such a way as to satisfy the boundary conditions (4.4) and 
(4.5) and the conditions at the origin and at «. We shall see that such solutions of 
(4.8) can be determined, which means that potential functions g(x, y) = Ref(z) satisfy- 
ing our conditions will be shown to exist. It will also be shown that the solutions ¢ 
of our problem behave on the water surface at ~ like C cos(x+a) in which C and a, 
the “amplitude” and “phase” of g, may have any values. Once the phase and ampli- 
tude at © are prescribed, however, the solution is uniquely determined. 

We proceed to solve the differential equation (4.8) and to fix the constants appro- 
priately. One integration can evidently be carried out at once to yield 


d i 
(i oe 1) f9 =—A—,; Areal. (4.9) 
dz z 


The additive constant which arises through the integration would be imaginary be- 
cause of the boundary condition (4.4); we have taken it to be zero since it would upon 
integrating (4.9) give rise only to an additive imaginary constant in the general solu- 
tion for f(z) and this in turn would contribute nothing to the real part ¢ of f(z). 
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A solution f;(z) of (4.9) for A =0 (i.e. a solution of the homogeneous equation) 
can be found which satisfies all of our conditions. The solution is 


fi(z) = Be-*, (4.10) 
with Ba real but otherwise arbitrary constant, as one can readily verify. The homo- 
geneous differential equation thus furnishes solutions of the problem which are 
bounded. The corresponding real potential function ¢g,= Ref;(z) is, evidently 

gi(x, vy) = Be’ cos x. (4. 10’) 
We observe that these solutions differ in amplitude but not in phase. They are, in 
fact, the non-singular solutions of our problem mentioned earlier in this section. 


Other types of solutions result from the non-homogeneous equation, i.e. for the 
case A ~0 and these will be singular at the origin. One solution of the non-homogene- 


~is gt 
f(z) = - Ae f — dt, (4.11) 
0 t 


in which the path of integration is taken along the positive real axis from + «, then 
along a circular arc about the origin, and then along a ray to the point 28, with 


ous equation is given by 








8, = —1, as indicated in Fig. 2. The integral evidently converges. However, it is 
y 
Z= a thy 
o -" _— 
= + co 
Z 
8, 


Fic. 2. Path of integration. 


necessary to add an appropriately chosen solution of the homogeneous equation in 
order to satisfy the boundary condition on the negative imaginary axis. (The bound- 
ary condition on the real axis is automatically satisfied since f(z) satisfies (4.9) with A 
real.) From (4.11) we find 


df —iz gt 1 
aoe Nae Ge a| - ies f —dt+—f, (4.12) 
dz us z 


and we are interested in the value of the expression on the right when 2 is a point on 
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the negative imaginary axis. One sees readily that the straight line portions of the 
path of integration (both of which lie on the real axis in this case) contribute purely 
real quantities to [;(e-‘/t)dt, but that the semicircular part yields an imaginary 
contribution of amount" plus ai. Thus the real part of the right hand side of (4.12) 
for z on the negative imaginary axis reduces to the real part of Awe~** there, since 1/z 
and ze~ ** are both pure imaginary on this axis and A is real. If, therefore, we add the 
function +Azie~** to the right hand side of (4.11) the result 


—iz e~ 
fez) = —A jew f — dt — rie*| (4.13) 
+2 t 


will be a solution of the non-homogeneous equation for which Redf,/dz=0 on the 
negative imaginary axis, i.e. it will be a function satisfying the boundary condition 
on the imaginary axis. 

Thus fi(z) and f2(z) as defined in (4.10) and (4.13) are two linearly independent 
analytic functions which satisfy the boundary conditions. We observe also that fi(z) 
and fo(z) behave in the prescribed manner at the origin: The function f(z) is regular 
there while fe(z) obviously has a logarithmic singularity. 

We have still to check the conditions at ». For f,(z) the conditions are obviously 
satisfied. To investigate the behavior of f2(z) as z—> we must consider the asymptotic 
behavior of the integral [&(e~‘/t)dt. We shall show in Appendix I that the integral 


possesses the following asymptotic representation :" 


“sf t : re 1 
— dt = 2a3 — te*| —+--- 
hee t Zz 


the dots representing higher order terms in 1/z. Assuming that such a development 
holds, then fe(z) clearly possesses the following asymptotic development : 


fo(z) S— Awie* = — Ane—*-*/?) (with A real). (4.14) 


The derivatives of fe(z) have essentially the same asymptotic behavior, as one can 
readily see. Thus all conditions at © are satisfied. The function ¢2(x, y) =Ref2(z) be- 
haves as follows at «: 
go(x, y) = — Ame” sin x. (4.14) 
The solutions ¢;(x, y) (cf. (4.10’)) and ge(x, y) are thus out of phase at © and we have 
therefore achieved one of our main objects. Our conjecture that solutions behaving 
like g, at © would result by imposing a logarithmic singularity at the origin proves to 
be correct. 
It is important to show that the set of all analytic functions f(z) which satisfy our 


boundary and regularity conditions is given by 
f(z) = filz) + folz) (4.15) 


with f;(z) and fe(z) the above defined solutions of (4.9) —/fi of the homogeneous equa- 
tion and fz of the non-homogeneous equation. This one proves as follows: Suppose 
that g(z) is any solution of our problem and set G(z) =f(s) —g(z), with f(z) given by 





entire value of the integral. 
4 The present case corresponds to the case 8, = —i of Appendix II, which in turn arises for k =2, n =2. 
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(4.15). It is clear that G(z) satisfies the same boundary conditions as f and g. Since f 
and g both satisfy (4.9) with the same value of the constant A it follows that 
(id/dz —-1)G=0 and the only solution of this equation satisfying our conditions is 
G=Ce-*, with Creal. Thus g(z) could differ from f(z) only by an additive real multiple 
of fi(s)—that is, g and f are really the same manifold of solutions. We observe that 
the set of solutions (4.15) contains two real constants A and B which are at our dis- 
posal. 

From (4.10’) and (4.14’) we conclude that real potential functions solving our 
problem and having any prescribed amplitude and phase at © can be obtained by 
superposition of g; and ¢g:, since A and B may be chosen arbitrarily, and since our 
boundary value problem for the function ¢ is linear and homogeneous. Conversely, 
it is evident that the constants A and B are uniquely determined when the phase and 
amplitude of any solution g at » are prescribed. Once this has been done the complex 
potential function f(z) is uniquely determined (cf. (4.15) and the remarks concerning 
it) and hence also the real potential function g(x, y). In other words, our solutions ¢ 
exist and are uniquely determined when we prescribe the phase and amplitude at ~. 

We now reintroduce the original space variables by replacing x and y in all rela- 
tions by mx and my, in which m satisfies the conditions o? = gm =27g/d (cf. Sec. 3) 
and @ is the circular frequency. At © our solutions ¢; and g2 have been shown to be- 


have as follows: 
gi = Cie™ cos mx, 


G2 = Cre™ sin mx, 
and consequently the quantity \ = 27/m is the “wave length” at ». This substantiates 
the remark made in Sec. 3 that the asymptotic relation between the wave length 
at © and the frequency is the same as that for water which is everywhere (i.e. for all 
values of x) infinite in depth. 

The standing wave solutions ¢; and ¢ will be discussed in detail in Sec. 6. 

5. Traveling waves over a sloping beach. In this section we shall generalize the 
method of the preceding section to yield solutions for waves on a beach which slopes 
at an angle 7/2n with the horizontal, with m any integer.’ The method we use is in 
principle exactly the same as for our previous case of a vertical cliff. The only differ- 
ences arise from the fact that the differential equation corresponding to (4.11) cannot 
be obtained quite so easily: it will be, in fact, a differential equation with constant 
coefficients of order 2n instead of one of order two. Naturally, the actual determina- 
tion of the desired solution will therefore become more and more complicated as n 
increases, i.e. as the inclination angle of the beach decreases. 

We formulate our problem at the outset in terms of the analytic function 
f(z) =e+. We seek such a function in the sector of angle 7/2n indicated in Fig. 3. 
The function f(z) should be regular in the interior of this domain, and have at most a 


18 The problem can be solved for other angles by similar methods, and probably also for any angle by 
extensions of the theory along the lines of the method of Sommerfeld used in diffraction problems. How- 
ever, it seems certain that such solutions would be very complicated and would involve functions not easily 
handled numerically with the tables of functions now available. The cases we discuss in this section are, 
it happens, amenable to numerical treatment. In any case, for angles less than 90°, it seems certain that 
the main features of the wave motion will be completely revealed through study of our special cases. For 
angles greater than 90°—that is, for overhanging cliffs or docks—new features could be expected to arise, 
and these cases deserve study. As we mentioned in the introduction, H. Lewy [6] has solved the problem 
for angles px/2n with p any odd integer and m any integer such that 2n >p. 
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logarithmic singularity at the origin, which we interpret to mean (cf. the remarks on 
this point in the preceding section) that |d*f(z)/ds*| <M;/|z|* for R=1, 2, - + - , 2n, 
with M, certain constants. At © we require that |®ef(z)| and |d*f(z)/dz*| for 
k=1, 2,--+, 2m should remain uniformly bounded when z— © in the sector.'® As 
a consequence, all derivatives of f(z) tend to zero along certain rays. On the boundary 
the conditions on g(x, y) =®ef(z) lead, as before, to the following conditions on f(z) 
at the boundary (cf. (4.4) and (4.5)): 








d 
ReL(D) - f(z) = te ( - po thesen? r )s = for s = rel) ¢ >0 (5.1) 
az 


MeL, (D) - f(z) 


d 

ae (i — i) J =(Q for z=2z>0. (3.2) 
dz 

By D we mean, of course, differentiation with respect to z. The boundary condition 

(5.1) should state that the derivative of g(x, y) normal to the bottom vanishes; that 

it does can be checked 








easily, for example by y 
inserting o+ for f, 
replacing d/dz by 0/0x, 
and using the Cauchy- 
Riemann equations. The 
condition (5.2) at the 
free surface is the same 
as (4.4). 

In the case treated 
in the preceding sec- 
tion, the operators L; 
and Ly» had the prop- - O Free Surface =X 
erty ReL,- Lof(s) =ReLe LLL Won ae 
-Iif(z) = 0 on both i Pag ting ia tad 
boundaries (cf. (4.6)). —— —— ——— —— 
This is, however, not Bottorh err 
the case for the corre- it ag 
sponding operators L, 
and Lo», defined in (5.1) 
and (5.2). It is neces- Fic. 3. 
sary for our purposes, 
in fact, to make use of a set of 2m linear operators Ly, Lo, - + +, Lena, Len, of which 
L, and Le, are the first and last members of the sequence, and which are so defined 
that Me(L,-Le-Ls- - - - Len-1-Lenf)=0 on both boundaries. It would be possible to 
derive these operators through geometrical constructions and arguments (involving 
essentially a succession of reflections in the lines re~***/?™, k=1, 2,---), but we 


prefer to give them at once and then examine them to see that they have the proper- 
ties we wish. They are defined as follows: 
{ (axD) if k is odd, 
aS ‘ ‘ 
Ua —1) if kis even, 


16 It would be possible to weaken this requirement considerably. 


(5.3) 
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in which the a, are the following complex numbers: 
ay = em ttlb/ntl)/2 k=1,2,---,2n. (5.4) 


It is useful to bear in mind the location of the points a, in the complex plane, as given 
in Fig. 4. These numbers lie on the unit circle spaced at equal angles 7/2n, all of 
them except a2, =7 having negative real parts. 

We show that these operators have the required properties. To begin with, we 
verify at once that the operators Z; and Le, as given by (5.3) are the same as those 
given in (5.1) and (5.2). We write: 


L(D)f = Ly-L2- ++ Lan: f(z) = (awD)(a2D — 1) + - (a2n—1D)(omD — 1) f(z). (5.5) 


Our object is to show that ReL(D)f=0 on both boundaries of the sector. We know 
that Re(a2,D —1)f(z) =0 on the real axis (condition (5.2)). We proceed to show that 
the operator P,(D) defined through (5.5) by L(D) =P,(D)-(as,D—1) has all of its 
coefficients real. It is clear that we may write the polynomial P,(D) as the product 
of two factors: P;(D) =P (D)P{' (D), with Pi (D) and P{’ (D) defined as follows: 


P{(D) = [are2n—1D?|[ascn—sD?] --- , 
P{’(D) = [(a2D — 1)(a2n 2D — 1)][(asD — 1)(an-sD — 1)] - - - 


That is, we separate the linear factors of P; into two groups, one containing all fac- 
tors for which & is even, the 


Y other all those for which k is 

bie ied odd; these two groups are 

Xn p én then arranged in the manner 

z indicated. From the defini- 

a a, tion (5.4) of the a, (cf. also 
J Zn Fig. 4) it is clear that a, 


=@2,-%, in which the bar 
over a quantity means that 
the complex conjugate of 
the quantity is taken, and 
also | ax| =1 for all k. Hence 
== X both ax-aen_x, and ap+aen_y 
are real numbers for all & 
={. 2. +++, 2a~—1, and 
hence all of the quadratic 
factors in P{ and PY’ obvi- 
ously have real coefficients. 
Since P;(D) contains an odd 
number of linear factors it 
follows that either P/ or P? 
will contain one unpaired 
linear factor, i.e. the factor 











Fic. 4. containing a,D. But since 
a, = —1 (cf. (5.4)) it follows 
that the polynomials P/ and P/’ have all their coefficients real and hence also the 
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polynomial operator P,(D). Consequently when z varies along the real axis we may 


write 


0 
ReL(D)f = Ps(—) SRe(anD — 1)f, 
< 


since D may be replaced by 0/dx in this case and the coefficients of the operator 
P,(0/dx) are all real. Therefore, in view of (5.2) we may write 


ReL(D)f = 0 on the real axis. (5.6) 


To prove that the same relation holds along the bottom, i.e. for z=re~*/*", we con- 
sider the operator P2(D) defined by L(D) =(a,D)P2(D). Along z=re-**/?" the opera- 
tor D=d/dz may evidently be replaced by the operator e**/"0/dr. By 0/0r is 
meant, of course, the directional derivative along the bottom line. If, then, D 
is replaced by e**/2"0/dr in P2(D) we observe that a2D, a3;3D,--++, a@2,D become 
00/Or, 020/Or, + + + , Qe»10/8r since ax41-e**/?" =a, from the definition (5.4) of the 
a,. It follows that the operator P:(D) along z=re‘*/?" may be replaced by an operator 
P;(0/0r) such that all coefficients of Ps; will be real. This can be seen in exactly the 
same manner as in the preceding case. We may write, therefore, along this line: 


0 0 
MeL(D) f(z) = RePs (: ) avyfee) = P; (—) weCaDyf(), 
r r 


0 
and, in view of the boundary condition (5.1), we see that 
ReL(D) f(z) =O on zs = re~**/?n, (5.7) 


Consequently, ReL(D)f(z) =0 on both boundaries. 
We now continue in the same fashion as in the preceding section by introducing 
the analytic function F(z) defined in our sector by 


F(z) = L(D)f(z), (5.8) 


and obtain a differential equation for f(z) by determining F(z) uniquely through use 
of the properties postulated for f(z). The properties of F(z) are the following: 1) Its 
real part vanishes on both boundaries of the sector, 2) It is either regular at the origin 
or it has a pole of order at most 2n, 3) | F(z)| remains uniformly bounded as z— 
in the sector. We know therefore that the derivatives of F(z) will tend to zero along 
any rays not parallel to the sides of the sector when |z|—>, by the Lemma intro- 
duced in Sec. 2. The first property has just been established, and the second and 
third follow immediately from our conditions on f(z) at the origin and at ~, through 
the definition (5.8) of F(z). 

Because of property 1) and the fact that we are working in a sector of angle r/2n, 
with » an integer, it is clear that F(z) can be continued over the boundaries of the 
sector by successive reflections to yield a single-valued function regular in the en- 
tire plane except possibly at the origin where it may have a pole of order 2 at most. 
In addition | F(s)| is bounded at » and the real part of F(z) vanishes on all rays 
through the origin for which z=re‘**/2", kR=1, 2, - - - , 4n. The only analytic function 
with these properties is F(z) =Ae,i/z2", with As, an arbitrary real constant (zero not 
excluded). (We rejected a possible additive constant in F(z) through use of the fact 
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that F(z) tends to zero along certain lines.) It follows therefore that the solutions f(z) 
of our problem satisfv the differential equation 


i 
L(D)f = (aD) (aD —_ 1)(a3D) (aD —_ 1) ie (@2n—1D)(ae,D — 1) f = Aan a ’ (5.9) 


with A», an arbitrary real constant. This differential equation may also be written 


in the form 


010305 * * * Aon-1D"(a2D — 1)(agD — 1)(agD — 1) -- - (@2,D — 1) f(z) = Aon — - (5.10) 

If we integrate (5.10) once with respect to z on both sides the only change 
on the left side is that D"* becomes D*—! while the right hand side becomes 
—Ao»,i/(2n—1)z*""-!+C, in which C, is an integration constant. But since | d*f/dz*| 
for k=1, 2, - - - , 2m tends to zero as | z| —c along any ray not parallel to the sides 
of the sector it follows that C;=0. The same process can be repeated—in particular, 
the successive integration constants will be zero for the same reason—until we ob- 


tain!” 





i 
(5.11) 


on 
~ 


(a2) — 1)(ayD — 1) - ++ (a2,D — 1)f = A 


with A real (but otherwise arbitrary) since the product aiasa5 - - - @2n-1 is real. The 
nth integration would also yield an additional constant on the right hand side of 
(5.11) which would not necessarily be zero, but which would be pure imaginary. This 
follows from the boundary condition Re(az,D —1)f=0 on the real axis, the fact that 
(a2D —1)(a4D—1) - - - (Qen-2D—1)=Pi’ (D) has only real coefficients (as we have 
seen a little earlier), and the obvious fact that ReAi/z”" (A real) is zero for real z. 
However, such an imaginary constant on the right hand side of (5.11) would clearly 
contribute to the general solution f(z) of (5.11) only an additive pure-imaginary con- 
stant which would contribute nothing to the real potential function g. Consequently 
we take this constant to be zero. 

Once we have the differential equation (5.11) the general procedure as well as the 
results are exactly analogous to those of the preceding Section 4: the arbitrary con- 
stants in the general solution of (5.11) can be fixed so that two types of solutions 
fi(z) and f2(z) are obtained which satisfy the boundary conditions'’® as well as the 
conditions at the origin and at «. The solutions f,(z) are obtained from the homo- 
geneous equation, while the f2(z) are solutions of the non-homogeneous equation. The 
solutions f;(z) are regular at the origin, while the solutions f2(z) have a logarithmic 
singularity there. All solutions f(z) of our problem are then given by f(z) =fi(z) +f2(z) 
and this set of solutions depends only on two real constants which occur as factors in 
fi and fe. At the behavior of f; and fe is such that the real potential functions 
i= Ref; and g2=Refe behave like Cye” cos (x-+a1) and Cy,e” cos (x+ a2) in which C; and 


17 The differential equation (5.11) was obtained by H. Lewy in a different way through reflecting in 
the bottom and then working in a sector of angle z/n instead of one of angle x/2n. At « Lewy assumes 
at the outset that the solutions behave like those in water of infinite depth, in contrast to the above 
treatment in which only the boundedness of certain derivatives is assumed. 

18 That this can be done is far from trivial, since our boundary conditions must be satisfied at all 


points on the straight lines composing the boundary of the sector. 
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C, are arbitrary, but a; and ag are fixed and differ by 7/2. It then follows that any solu- 
tion ¢ of our problem is uniquely determined once the phase and amplitude of g at 
are prescribed. 

The general solution of (5.11) containing m arbitrary constants is of course ob- 
tained by straightforward and elementary methods. However, the determination of 
these constants in order to satisfy the boundary conditions is not entirely trivial, 
particularly in the case of the solution fe(z) of the non-homogeneous equation. In 
Appendix I we discuss the method of determining the integration constants in such a 
way as to satisfy the boundary conditions; in the present section we simply give 
the results. The solutions are also seen to satisfy the conditions at the origin and at . 

In the homogeneous case, the solution is 


w n 
fi(z) = = Dice, km i,2,---,m. (5.12) 
‘ (n—1)!/n pea 


The numerical factor before the summation sign is chosen for convenience in later 
computations. The constants c;, and 8; are the following complex numbers: 


Bi = e'™ k/n+1/2) (5.13) 
neee T 2r (k — 1) 
Cc. = ett (mt) /4-k/2] Cot — cot —--- cot ———— for k= iH ~% inert. 2 
2n 2n 2n (5.14) 
C1 = Cy 


By comparison with (5.4) we note that 8, =a. (The bar on é, and & means that the 
complex conjugate of these quantities is taken.) The constants c, are uniquely de- 
termined (see Appendix I) within a real multiplying factor. 

As |z|—>© in the sector, all terms clearly die out exponentially except the term 
for k=n, which is c,e~*, since all 8,’s except 8, have negative real parts. Even the 
term for k=n dies out exponentially except along lines parallel to the real axis. (The 
value of c,, by the way, is e~**‘"—)/4 since the cotangents in (5.14) cancel each other 
for k=n.) This term thus yields the asymptotic behavior of f;(z) : 

Tv 
fiz) —— =e, (5.15) 
; (n — 1)!V/n 


The general solution f2(z) of the non-homogeneous equation (5.11) when the real 
constant A is set equal to one is as follows: 


n zBk e-? 
f(z) = Saal em f — di — ries, (5.16) 
k=1 + t 


The §;’s are defined by (5.13); and the a,’s are defined by 
an = cr/(n — 1)!V/n, (5.17) 


that is, they aré a fixed multiple (for given m) of the c,’s defined by (5.14). The a,’s, 
like the c,’s, are uniquely determined within a real multiplying factor. The path of 
integration for all integrals in (5.16) is indicated in Fig. 2. That the points 28; lie in 
the left half of the complex plane (as indicated) can be seen from our definition (5.13) 
of the 6; and the fact that z is restricted to the sector 0 Sarg zSa/2n. 

The behavior of f2(z) at «© of course depends on the behavior of the integrals in 
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(5.16). In Appendix II it is shown that these integrals behave asymptotically as 


follows: 


e~ Pk / 1 T 
_ ( bee), — < arg 28; Sm, 
Bk ‘di B. \z 2 = 
~ (5.18) 
oY iv t é 2Bh 1 3 
2ri — ( +--+: ). ew < arg 26:3 
By Zz 2 


Once this fact is established it is clear from (5.17) and (5.16) that fo(z) behaves asymp- 
totically as follows: 
T 

f(s) ~- — iad (5.19) 

: (n — 1)!V/n 
since the term for k=n dominates all others (cf. (5.18)) and arg 28, >~7 in this case. 
Comparison of (5.19) with (5.15) shows that the real parts of fi(s) and fo(z) would be 
90° out of phase at «. 

That the derivatives of fo(z) behave asymptotically in the same fashion as f.(z) 
itself is easily seen, since the only terms in the derivatives of a type different from 
those in (5.16) are of the form 6,/z*, k an integer 21. Finally, it is clear that fo(z) 
has a logarithmic singularity at the origin. Hence f;(z) and f2(z) satisfy all require- 
ments. Just as in the 90° case (cf. the preceding section) it can now be readily seen 
that f(s) =dif;(z) +bcfe(s), with b; and be any real constants, contains a// standing wave 
solutions of our problem. The proof that the real potential function g=®ef(z) is 
uniquely determined once the phase and amplitude at ~ are fixed can also be carried 
out exactly as in the previous section for the 90° case. 

The relations (5.15) and (5.19) yield for the asymptotic behavior of the real po- 


tential functions ¢g, and ¢2 the relations: 


vs n—1 
gi(x, y) = Refi~ e” cos (« + — (5.20) 
: (n — 1)!/n 4 
T : a1} 
g2(x, y) = Refe~ ev sin { x + —— (5.21) 
(n — 1I)lV/n 4 


when it is observed that c, =e~'*‘"~"/*, It is now possible to construct either standing 
wave or progressing wave solutions which behave at ~ like the known solutions for 
steady progressing waves in water which is everywhere infinite in depth. In particular 
we observe that it makes sense to speak of the wave length at © in our cases and that 
the relation between wave length and frequency satisfies asymptotically the relation 
which holds everywhere in water of infinite depth. For this, it is only necessary to 
reintroduce the original space variables by replacing x and y by mx and my, with 
m=o"/g (cf. Sec. 3), and to take note of (5.20) and (5.21). 
Finally, we write down a solution P(x, y;t) which behavesat @ likee”cos (x+t+a), 


i.e. like a steady progressing wave moving toward shore: 
#(x, y; 2) = A[gi(x, y) cos (t + a@) — ¢g2(x, y) sin (t+ a)]. (5.22) 


The solutions (5.22) will be discussed numerically in Sec. 8 for the case of a beach 


sloping at an angle of 6° (i.e. for the case n=15). 











| 
| 
| 
| 
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6. Numerical discussion of standing wave solutions for 90°, 45°, and 6° slopes. 
In this section we give graphs of the two types of standing wave solutions for the case 
of a vertical cliff (cf. Sec. 4) and for bottom slopes of 45° and 6°. We continue to 
make use of the dimensionless variables of Sec. 3. In particular, it should be re- 
called that the variable x means the distance from shore divided by 4/27, in which A 
is the wave length at infinity. In other words the quantity x is proportional to the 

wave length at ~. 

In the case of the vertical cliff, or 90° case, two standing wave solutions are 

given by 


| 
(x, y; 4) = mwe*'e” cos x (6.1) 


=cos £ * sin & 3 
b.(x, y; 2) e* tev | cos rf - dé + sin vf "— dé + m sin «|, (6.2) 
wo E i) 


readily verify, either from (4.10) in Sec. 4 or from (5.12) in Sec. 5 with 


and 


As one can 
n= 1. 
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Fic. 5. Standing waves for a vertical cliff. 
#,(x, 0,0) —-—- x =distance from shore/(\/27) 
@,(x, 0, 0) - \=wavelength at © 
7 sin x ----- 


The functions ®,; and ®, for y=0 and t=0 are plotted"® in Fig. 5 together with the 
function x sin x, which yields the asymptotic behavior of ®: on the surface. The most 


19 The curves for &; and 4, differ from the corresponding surface elevations m and m (cf. (1.4)) only 


by a phase shift and a constant multiplier. 
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notable feature of the curves of Fig. 5 is that the standing wave solution ®2, although 
it tends to become negative infinite as x0, nevertheless differs little from the func- 
tion 7 sin x until a point within a distance from the cliff of a fraction of the wave 
length at © is reached. In other words, the asymptotic development of ® for x large 
holds with good accuracy until x is rather small. In addition, the last maximum of ®, 
(i.e. the crest nearest to the cliff) has an amplitude slightly Jess by a little more than 
10 per cent than the amplitude at ©. The distance between the two zeros of ®, 
nearest to shore is also 10 per cent less than half the wave length at ~. 

An interesting additional fact which is not difficult to prove is that the velocity 
of the water does not decrease exponentially with the depth y (as it does when no cliff 
is present), but only like 1/y. The presence of the cliff thus reinforces the velocity. 
(The velocity of the water along the cliff is given by Me(¢df(z)/dz) with f(z) defined 
by (4.15); the result indicated follows by calculating idf/dz and using the known 
asymptotic behavior of the complex integral which occurs.) 

The standing wave solution #,, ,. for the case of a beach sloping at 45° are ob- 
tained from (5.12) and (5.16) with »=2. For ®, we write 


Tv : 
®,(x, vy; 4 = = eRe le* te-2 + e-inlte—iz| (6.3) 
V2 


The unbounded standing wave solution is given by 


et! f ; 2 gre ; 
@o(x, y; f) = —= Re ye*r!4 otf —— dt — ric” 
sz e? 
+ e'*/4 [ «f — dt — rie “lt . (6.4) 
p t 


The surface values of ®; and ®2 for t=0 are plotted in Fig. 6. These curves are ob- 
tained by using tables of the functions C;, S;, and E,*° computed by the Mathematical 
Tables Project [7]. In fact, , can be written in the form 

e*tey 
$.(x, y; 4) = — = Ca) [sin x — cos x] 

V/2 


T ; \ 
_ E + S;(x) [cos x + sin «| —€ lied | ; (6. 4’) 


At «, (x, 0; 0) behaves like (#/+/2)(cos x—sin x) while ®:(x, 0; 0) behaves like 
—(r, ‘~/2 \(cos x+sin x). 

The same general comments can be made about the curves of Fig. 6 as were made 
for those of Fig. 5. In particular, the minimum of ®, at x =i is about 10 per cent less 
in numerical value than the amplitude of ®, at ~, but differs very little from its 
asymptotic representation until x is quite small. The regular solution ®, attains a 
maximum on shore which is+/2 times the amplitude at . The distance between suc- 
cessive zeros shortens on approaching shore, as in the preceding case, but the shorten- 
ing is now more pronounced. 





20 We use i unconventionally as a subscript to avoid confusion with i= 4/ —1; in M.T.P. [7] the nota- 


tion Ci, Si, and Fi is used. 
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Fic. 6. Standing waves for a 45° bottom slope. 
,(x, 0, 0) —-—- x =distance from shore/(\/27) 


(x, 0, 0) — \=wavelength at 
—r/+/2(cos x+sin x) ---- 


Finally, we describe the two types of standing wave solutions for the case of a bot- 
tom slope of 6°. These solutions are obtained from (5.12) and (5.16) by taking »=15. 
The regular standing wave ®, is taken in the form (cf. (5.12)): 


®i(x, y;) = e* Refi (z) (6.5) 


with f;(z) defined by 


15 
fi(z) = Do cree, 
k=1 
By = eft 1541/2), 
T (k—1)x 
t a 


Cc. = e7itk/2 cot —--- co 


30 30 


(6.6) 
» kA, 


Oo=—4, 


The fifteen quantities c, are alternately real and pure imaginary and they vary in 
absolute value from 1 to approximately 819. For large values of x the function 
®,(x, 0; 0) behaves as follows (cf. (5.20)): 


@,(x, 0; 0) ~ x sin x. (6.7) 
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The standing wave solution #2, which is infinite at the shore line, is defined by 
@(x, yj 1) = — 141/15 e* MRefo(z) (6.8) 


with fe(z) defined by 
1 > f zBk e t 
fo(z) —————— Ck com f - dt — ries (6.9) 
. imo tt s..t 
in which the c;, and f; are defined as in (6.6).- The path of integration is shown in Fig. 2 
of a preceding section. This solution behaves for large x and for y=0, t=0 as follows 
(cf. (5.21)): 
$2(x, 0; 0) ~ wr cos x. (6. 10) 


In order to determine ®, numerically the function 


oo e t 
E(z) = ef . dt 
P l 


was computed for values of z on the rays z=re**(*/5+"2), kR=1,2,---,7. 
The function E(z) defined by (6.11) has been tabulated for values of z in the second 
quadrant by the Mathematical Tables Project [8]. However, the entries in this table 
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Fic. 7. Standing waves for a 6° bottom slope. 
(x, 0,0) —-—- x =distance from shore/(A/27) 
@,(x, 0, 0) \=wavelength at « 
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are for a rather wide rectangular net of points in the z-plane. We therefore found it 
necessary to compute E(z) along the above described rays as follows: For r= | s| = 3.00 
and r=3.50, we interpolated in the M.T.P. tables of E(z) by means of a Taylor series 
expansion about the nearest tabulated points. (The derivatives of E(z) are easily 
calculated.) The values at r=3.00 and 3.50 were checked by using them to compute 
the values at r=3.25. Power series developments were thus obtained for points along 
each of the seven rays for the range 1.00 Sr $7.00 at an interval of 0.25, which means 
that more than .100 power series were derived. A table was then constructed by using 
these series to obtain values at intermediate points. Since integrals of this form would 
always occur in solving differential equations having constant coefficients and rational 
right hand sides, it seems worth while to include this table as an Appendix. The table 
is believed to be accurate within one unit in the fourth decimal place. The table 
would, of course, also be useful if computations were to be made for beach slopes at 
angles other than 6°. 

Fig. 7 shows the surface values of ©, and ®, near shore, while Figs. 8 and 9 are 
graphs of ®, and ® together with their asymptotic representation for large x. Again 
the same general comments are in order as in the preceding cases, except that it is 
now necessary to go further out from shore to obtain close agreement with the 
asymptotic solution in deep water. This is not very surprising since the depth of the 
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Fic. 8. Comparison of standing wave solution for a 6° slope, #,, with 
its asymptotic limit sin x (see Fig. 7). 
(x, 0,0) —-—- x =distance from shore/(A/27) 
msin (x) - — \=wavelength at 
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Fic. 9. Comparison of standing wave solution for a 6° slope, 42, with 





its asymptotic limit # cos x (see Fig. 7). 
,(x, 0,0) ——— x =distance from shore/(A/27) 
wcos (x) ---- \=wavelength at 


water at x=7 (the largest x for which we have numerical values) because of the small 
slope is less than 1/8 of the wave length at «. Comparison of Fig. 7 with Figs. 5 
and 6 shows that the distance between successive zeros near shore has now shortened 
very materially as compared with the preceding cases. In fact, this effect would be- 
come more and more pronounced with decrease in the slope.*! We observe also that 
the relative maximum of ®, at x =5.4 and the relative maximum of ®,; at x=6.5 are 
less than the amplitude at ©. Progressing waves for the case of a 6° slope will be 
discussed in Sec. 8. 

7. Shallow water theory. In some gravity wave problems it is possible to obtain 
an accurate approximation to the exact linear theory by relatively simple means. 
Such an approximate theory, which is commonly referred to as the linear shallow 
water theory, has long been the basis for the theoretical study of the tides in the 
oceans and large estuaries. In this section we give a derivation of the shallow water 
theory for water of variable depth; in the next section we shall compare the results 
from both theories numerically in the case of progressing waves in water of constant 
depth and over a beach with a 6° slope. In what follows we consider only the two- 
dimensional case. We revert also to the original space and time variables. 


21 See the next section, where an approximate theory valid for water in which the depth is small 


compared with the wave length is discussed. 
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The exact linear theory requires the determination of a potential function 
®(x, y; t) in the region indicated in Fig. 10). As before, y=0 is the original undis- 
turbed surface of the water and y= —h(x) is the bottom profile. The elevation of the 


y 


| 














Bottom 


Fic. 10, 


free surface in the course of the motion is denoted by n(x, t). The conditions to be 
satisfied by ® are (cf. Sec. 1). 


ft t= 8 for 02 y= — h(x) (7.1) 

Pi, = — go, for y=0 (7.2) 

&,= —h,?, for y= — h(x). (7.3) 

For the purposes of the present section it is not necessary to formulate conditions 


at «. The surface elevation 7 is given by 
bd , 
n(x,t)=—-®| . (7.4) 
g y=0 
In the following derivation it is convenient to denote quantities evaluated for y =0 
by a bar over the quantity, and for y= —h(x) (i.e. at the bottom) by a bar under the 
quantity. Thus conditions (7.2) and (7.3) could be written in the form $= —g@y 
and ®,= —h,®,. 
We may write 


0 
f &,dy = $, — o, = %,+ h.%, (7.5) 


A(x) 
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from (7.3). On the other hand we may also write 
0 0 0 0 
f $,,dy = - f ?, zdy ese - —{ ?dy + h,?,; (7.6) 
—h(z) —h(2) Ox J —n(2) 


by using (7.1) and observing that the lower limit h(x) depends on x. By equating the 
right hand sides of (7.5) and (7.6) we find 


ra] 0 
&, = —-- f ?.dy. (7.7) 
Ox —h( 2x) 
We consider next the relation 


0 0 
f &.dy = he, -f yPidy (7.8) 
a h 


h(x) 


~ 


obtained through integration by parts, and also the relation 


0 
f h(x) Pzydy = h(x), — h(x) %,. (7.9) 


h(x) 


The quantity ®, in (7.8) can be eliminated by using (7.9), and this in turn lead- 


through use of (7.7) to 


fa) 0 
?, waa | maya, -f (y + i) Paty] 
Ox —hA(z) 


— (h®,), +f (vy + h) ®rey + hePy|dy. 
—h 


From the boundary condition (7.2) we may replace 6, by —1/gé,; to obtain 
Be = 0 
®,, = (gh®,). — ef [(y + h) Pry + hePry|dy. (7.10) 
-h 


Up to now we have made no assumptions in addition to those made in deriving the 
exact linear theory—we have simply carried out integrations with respect to y in 
such a way as to yield (7.10). We now assume that ®,,, and ®,, are bounded for all x 
and t, and for —h<Sy0, and that h, is small of the same order of magnitude as h. It 
follows at once that (7.10) may be written in the form 


®,, = (gh®,), + O(h’) (7.11) 


in which O(h?) is a quantity of order h?. Thus, if the depth / and slope h, are suffi- 
ciently small, the surface value @ of the potential function (x, y; t) should be given 
with good approximation by the differential equation 


®,, = (gh®,);. (7.12) 


In the case h=const. (7.12) is the wave equation in one space variable with 
c=(gh)'/? as propagation speed. Equation (7.12) is the differential equation of the 
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linear shallow water theory.** The above derivation of the theory is due to F. John. 
Instead of ® we may introduce the surface elevation n(x, 1) =@,/g (cf. (7.4)) as 
dependent variable; this leads easily to the equation 
nee = (ghnz)z. (7.13) 
It is possible to derive the shallow water theory in such a way as to obtain in 
place of (7.10) a relation in which the integral on the right hand side is replaced by a 
power series in h which converges if the initial conditions (i.e. conditions for t=0) 
are assumed to be sufficiently smooth. That is, it would not be necessary in this ver- 
sion of the theory to assume at the outset that ®,,, and ®,, are bounded for all time. 
We proceed to derive standing wave solutions of (7.12) for the case in which 


h(x) = qx, g = const., (7.14) 


i.e. for the case of a uniformly sloping bottom profile.** For this purpose it is con- 


venient to set 


@(x, t) = ef @'t+9Z (2), (7.15) 
The function Z(x) satisfies 
re m o 44 
(aZ2)2 + —Z = 0, m=—) (7.16) 
g g 


in which the quantity m has the same definition in terms of o as in the exact theory of 
the preceding sections (cf. Sec. 3). 
The differential equation (7.16) is a Bessel equation with the general solution 


mx mx 
Z(x) = AJo (24/™*) + By, (24/™*), (7.17) 
q q 


Jo and Yo are the regular and the singular Bessel functions of order zero respectively ; 
thus Y, has a logarithmic singularity at x =0. For large values of x these functions are 


well known to behave as follows: 


2 The derivation often given for this theory (cf. for example [4] p. 254) starts with the assumption 
that the pressure p is determined by the same relation as in hydrostatics, i.e. p=g(n—y). One would be 
inclined to think that such a relation would be on the whole more accurate the deeper the water since all 
motions die out in the depth, but it is easily seen on the basis of the simplest examples that the approxi- 
mate theory cannot be accurate in sufficiently deep water: The exact solutions for steady periodic waves 
in water of constant depth yield waves in which the wave length depends essentially on the frequency, 
but the steady waves given by (7.12) in the case h=const. are without dispersion. In other words, the 
derivation of the approximate theory by means of the hydrostatic assumption is open to criticism since 
it does not indicate clearly the essential role played by the depth in determining the accuracy of the ap- 
proximation. Cisotti has given another derivation of the shallow water theory (cf. [10] pp. 357, 379) which 
does not start with the hydrostatic assumption. 

23 This theory is, of course, well known. See, for example [4], p. 291. 
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(7.18) 


With this theory it is therefore not possible to obtain either standing waves or pro- 
gressing waves with non-vanishing amplitudes at ©, in sharp contrast to solutions 
given by the exact theory. Also, for large values of x the wave length (defined as the 
distance between successive nodes, say) as given by the approximate theory would 
be roughly 27(gx/m)*/* and would therefore increase indefinitely with x. This is also 
in marked contrast with our exact solutions, in which the wave length at © tends to 
a constant. 

Nevertheless, it is possible to reintroduce the time factor and obtain for n(x, ¢) 
solutions which have the form of progressing waves traveling toward shore, as follows: 


@(x, t) = A[cos (ot — €)¥o(2\/mx/q) + sin (ot — €)Jo(2y ‘mx/q) |. (7.19) 


One could expect that such a solution might furnish in some cases at least a fair ap- 
proximation to the exact solution for a not too large range of values of x. In particu- 
lar, we note that the singularity at the origin (i.e. on the shore line) is of the same type 
as in the exact theory. In the next section a numerical comparison with the exact 
theory will be made. 

8. Comparison of exact and shallow water theories. In the present section we com- 
pare the results obtained using the shallow water theory derived in the preceding 
section with those of the exact theory. We consider progressing wave solutions first 
for the case of water of uniform but finite depth and then for the case of a bottom 
slope of 6°. In this section it is preferable to use the original independent variables 
rather than the dimensionless variables of Sec. 3. 

In the case of water of uniform depth / the velocity potential P(x, y; ¢) is given 
by (cf. [4], pp. 363-368) 
cosh m(y + h) 

&(x, y; 1) = A ——————— sin (mx + ot), (8.1) 
cosh mh 

in which h represents the constant depth of the water. The undisturbed surface of the 
water is the line y=0. This potential function evidently satisfies the condition 
d@/dy=0 at the bottom y= —h. The free surface condition (1.2) is also satisfied if 


o? = mg tanh mh, (8.2) 


or as we may also write 





m?*h? 


“ 


If we introduce the phase velocity c=a/m we may write 


al 
a gh(1-7—+---) (8.4) 
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Val Q)+ | @. 
ee eee aS 3 
. a% & , 
if the wave length \=27/m is used. The relations (8.2) to (8.5) characterize the type 
of dispersion encountered in this case.*4 If the ratio of wave length A to depth h is 
large, the relation (8.5) for the phase velocity c becomes 

c= vVgh. (8.6) 


In fact, if \/h=5, c=.82(gh)'/?; while if \/h=10, c=.94(gh)"/*. 
In the same case of constant depth the shallow water theory gives the following 
approximation (x, ¢) to the surface value ®(x, 0; ¢): 


or 





@ =A sin (mx + ot), (8.7) 
$(x0;t) 
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Fic. 11. Progressing waves for a 6° bottom slope. (Exact theory.) 
A\=wavelength at 
x =distance from shore/(A/27) 
This is a solution of (7.12) only if 
o?/gh = m?, (8.8) 
i.e. if the phase velocity c=oa/m is given by 
c= VV gh. (8.9) 


* This theory appears to be due to Airy [1], who published it in 1845. 
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In other words, the approximate solution (8.7) furnished by the shallow water 
theory yields the phase velocity with good accuracy only if the ratio of the wave 
length to the depth is sufficiently great, as we see by comparing (8.9) with (8.5). 
In fact, (8.9) yields ¢ correct, within about 5 per cent, only if the wave length is more 
than ten times the depth. For steady waves, therefore, the approximate theory is ac- 
curate only if the water is shallow compared with the wave length. For this reason the 
approximate theory is sometimes referred to as the long wave-shallow water theory. 

We turn next to a comparison of the results of the two theories for progressing 
waves over a beach sloping at 6°. The progressing wave solution of the exact theory 
is given by (5.22), which behaves at © like 7 sin (mx+ot). Graphs of the numerical 
solutions for times of =0, 7/4, and 7/2 are shown in Fig. 11. (Again we note that 
the dimensionless variables of Section 3 are not used.) We observe that the “ampli- 


tude” of the progressing wave increases as the wave moves from the point 27x/A=6 


$(7,t) 
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Fic. 12. Progressive waves for a 6° bottom slope (shallow water theory). 
\= wavelength at © in exact linear case 
x = distance from shore/(A/27) 
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toward shore. However, the maximum at this point is 6 per cent less* than the am- 
plitude at ~. 

In the preceding section it was already explained that the shallow water theory 
could not, in principle, furnish a good approximation to the progressing waves in the 
present case over the whole range of values of x since the amplitude tends to zero at © 
in this theory. We have therefore chosen to make a comparison by assuming that the 
circular frequency ¢ is the same for both theories and that both solutions have the 
same relative maximum for ot=2/4 and 2rx/A=6. This yields for the approximate 
solution (x, t) the following (cf. (7.19)): 


mx’ 
P(x, i) = — 14.5| cos (at — .98)Vo (24/7) 
— q 
mx 
+ sin (ot — .98)Jo (24/™*)]. (8.10) 
q 
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Fic. 13. Comparison of surface values from both theories. 
#= Exact linear theory \=wavelength at © (exact theory) 
@=Shallow water theory ———— x=distance from shore/(\/27) 
Phase and amplitude are the same at (x=6, ot = 2/4) 
Frequencies are equal 








% The author is told that this effect has been observed experimentally. 
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Fig. 12 is a graph of @ for ot =0, 7/4, #/2. The curves have the same general appear- 
ance as those of Fig. 11, but the amplitude increases somewhat more rapidly in the 
present case. This is brought out more directly by Fig. 13, which gives graphs of ® 
and @ for ot =1/4. The two theories agree fairly well as the shore is approached, al- 
though the amplitude given by the approximate theory at 27x/A~2 is about 15 per 
cent greater than that furnished by the exact theory. We have at our disposal the 
information necessary to describe how the continuation of the curves of Fig. 13 would 
appear for values of 27x/X greater than 7, since both curves are given with good ac- | 
curacy in this range by their respective asymptotic representations. Thus ® will be 
very closely the same as @ sin (27x/A-+7/4)—that is, the exact solution would be | 
one having an almost constant amplitude equal to t—while the amplitude of @ de- 
creases like 1/(27x/d)~4 (cf. (7.19)). If, therefore, we were to compare the amplitudes 
at 2xrx/h=24, the amplitude of @ would be about 40 per cent less than that of ®, 
while at 27x/ = 12 the error would be about 20 per cent. On the other hand, it should 
be stated that at 2rx/A=7 the depth of the water for the 6° slope is somewhat less 
than 1/8 of the wave length. Thus the shallow water theory might be expected to 
yield fairly accurate results from this point in toward shore, but it could not be ex- 
pected to do so for 2rx/A=12, much less for 27x/d = 24, on the basis of our discussion 
above for the case of water of uniform depth. 

It is of some interest to compare the phase velocity c furnished by the two theories 
for 2rx/X<7. This was done by calculating the position of the zeros and the maxima 
and minima for a series of closely-spaced values of ¢; the velocities were then obtained 
by taking difference quotients. The results of such calculations for 27x/AS7 are 
shown in Fig. 14. The asymptotic value for c as furnished by the exact theory is indi- 
cated. Up to a distance of about a wave length from shore the two theories yield the 
same phase velocity, but from this point on the phase velocity predicted by the 
shallow water theory is too high and becomes more and more inaccurate as the dis- 
tance from shore (and therefore also the depth of the water) increases. At a distance 
of about 3 wave lengths from shore the phase velocity as given by the shallow water 
theory is in error by 10 per cent. 

We remark, finally, that the curves of Fig. 14 can not be distinguished (to the 
scale used there) from the curves 





_ gr 2rh 
c=vVgh, and c= — tanh — (8.11) 
2r nN 


in which h is the depth corresponding to any given x and ) is the wave length at «. 
In other words, the phase velocities actually computed by using the two theories for 
2xx/d S7 are practically identical in each of the two theories with those which would 
be obtained by calculating c for each value of x through use of the corresponding 
theory applied to water of the correct uniform depth. 

9. A problem in three-dimensional wave motion. Except for the present section, 
we consider in this paper only motions which are the same in all planes parallel to an 
x —y-plane, and which therefore can be treated by using functions of a single complex 
variable. (The exact linear theory is, of course, in question here.) For motions which 
depend essentially upon three space variables it is not possible to make use of com- 
plex functions, but it is possible to extend the basic idea of the method used in the 
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two-dimensional cases to surface wave problems in three dimensions. In this section 
we illustrate the method by treating the problem of progressing waves in an infinite 
ocean bounded on one side by a vertical cliff—in other words, the same problem as 
that of Sec. 4 except that we no longer require the waves to move with their crests 


parallel to the shore line. 


on 


27 


allow Water Theory 


a 6 8 /0 /2 ot / 





Fic. 14. Phase velocity for progressing waves over a 6° sloping beach 
x=distance from shore/(A/27) 
\=wavelength at © (exact theory) 
c=phase velocity. 


We seek solutions P(x, y, 2; ¢) of Viz.y.2.P=0 in the region x 20, yS0, —-~7 <z<@ 
with the y-axis taken normal to the undisturbed free surface of the water and the 
z-axis”® taken along the “shore,” i.e. at the water line on the vertical cliff x=0. 
Progressing waves moving toward shore are to be found such that the wave crests 
(or other curves of constant phase) at large distances from shore tend to astraight 
line which makes an arbitrary angle with the shore line. For this purpose we seek solu- 
tions of the form 


B(x, y, 2; t) = efltthetPdo(x, y) (9.1) 

%* It has already been pointed out that functions of a complex variable are not used in this section, 

so that the reintroduction of the letter s to represent a space coordinate should cause no confusion with 
the use of z as a complex variable in earlier sections 
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that is, solutions in which periodic factors in both z and ¢ are split off. Of course, the 
function g(x, y) is not, as in our previous cases, a potential function; instead, it satis- 


fies the differential equation 


: ; (9.2) 


—-—g=0 for y=0. (9.3) 


which implies that the dimensionless space and time variables of Sec. 3 (including now 
z as well as x and y) are assumed at the outset. The condition at the cliff, is, of course, 
dg 
— 0 x = 0. (9.4) 


a for 


Ox 
At the origin x =0, y=0 (i.e. at the shore line on the cliff) we require, as in former 
cases, that g should be of the form 


g=glogrt+y, rX1, ' (9.5) 


for sufficiently small values of r = (x?+ y?)"/*, with g and ¢g certain bounded functions 
with bounded first and second derivatives in a neighborhood of the origin. The func- 
tions @ and ¢ should be considered at present as certain given functions; later on, 
they will be chosen specifically. 

For large values of r we wish to have ®(x, y, z, t) behave like eve *('+*#+2+8) with 
k?+a?=1 but & and a otherwise arbitrary constants, so that progressing waves tend- 
ing to an arbitrary plane wave at © can be obtained. This requires that g(x, y) 
should behave at @ like e¥e*(°7+#2) because of (9.1). However, it is no more necessary 
here than it was in our former cases to require that g should behave in this specific 
way at ~; it suffices in fact to require that 


Jolt] ee|+leu| <M for r> Re (9.6 


| 


i.e. that g and the two derivatives of g occurring in (9.6) should be uniformly bounded 
at ©. As we shall see, this requirement leads to solutions of the desired type. 

We proceed to solve the boundary value problem formulated in equations (9.2) 
to (9.6). The procedure we follow is analogous to that used in the former cases in every 
respect. To begin with, we observe that 

0/90 

— (" ~ i)e = 0 forbotha=0 and y=0, (9.7) 

Ox \Oy : 
because of the special form of the linear operator on the left hand side together with 
the fact that (9.3) and (9.4) are to be satisfied. A function (x, y) is introduced by the 


--(- 1 9.8 
ahs a i 


relation 
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The essential point of our method is that the function is determined uniquely within 
an arbitrary factor if our function ¢, having the properties postulated, exists. Further- 
more, ¥ can then be given explicitly without difficulty. The properties of w are as 
follows. 
1) w satisfies the same differential equation as ¢, i.e. equation (9.2), as one sees 
from the definition (9.8) of y. 
2) Wis regular in the quadrant x >0, y <0 and vanishes, in view of (9.7), on x =0, 
y <0 and y=0, x>0. Hence y can be continued over the boundaries by the 
reflection to yield a continuous and single valued function having continuous 
second derivatives ~.z and W,, (as one can readily see since Vy—k*y=0, and 
¥ =0 on the boundaries) in the entire x, y-plane with the exception of the origin. 
(Here we use the fact that our domain is a sector of angle 1/2.) 
At the origin, y has a possible singularity which is of the form (x, y)/r?, with 


3) 
y regular, as one can see from (9.5) and (9.8). This statement clearly holds for 
the function y when it has been extended by reflection to a full neighbor- 
hood of the origin. 

4) The condition (9.6) on ¢ clearly yields for y the condition that is uniformly 


bounded at © after y has been extended to the whole plane. 


Thus v is a solution of V*y —k*y =0 in the entire plane which is uniformly bounded 
at o. At the origin Y=yY/r? with J a certain given regular function (J=0 not ex- 
cluded). In addition, Y=0 on the entire x and y axes. It is not difficult to prove that 
the solution of this problem is unique.?? 

On the other hand, a solution y of the problem for a special function®* J charac- 
terizing the singularity at the origin is readily given in polar coordinates (r, 8); it is 


W(x, y) = AiH; (ikr) sin 20, r=/e8+y, OSkS1, (9.9) 


in which H$” is the Hankel function of order two which tends to zero as r— ©, and A 
is any real constant. The function wy has real values when r is real. (The notation given 
in Jahnke-Emde, Tables of Functions, is used.) One can readily verify that this func- 
tion really does satisfy all conditions imposed on y. For our purposes it is of advantage 
to write the solution y in the following form: 


27 One way to do sois the following: The difference, V, of two solutions would have all of the proper- 
ties of y except that it would be regular at the origin. If Vo is the value of ¥ at any point (xo, yo) in the 
plane, then it is known (see, for example, Courant, Hilbert, Methoden d. Math. Phys., Bd. II, S.261) that 
the mean value formula 

Wo: Jo(ikR) = M 
holds. Here M is the mean value of W over any circle of radius R and center at (xo, yo). The function Jo 
is the regular Bessel function of order zero. If Ris chosen large enough M remains less than a certain con- 
stant since WV is uniformly bounded at «. On the other hand, Jo(kR) behaves for large R like e*®(24kR)~? 
(see Jahnke-Emde, Tables of Functions, p. 138) and hence as R-+, Wo would tend to zero. But since 
Vois independent of Rit follows that Vo is zero at any arbitrary point (xo, yo). Thus ¥ =0, and the unique- 
ness of the function y is proved. 

28 Our uniqueness theorem is less general in the present case than in the earlier cases since we pre- 
scribe the singularity at the origin so specifically in the present case. 
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2 
y= Pre Hy (ikr), r=V/a®+y?, (9.10) 
Oxdy 
in which A is any real constant and H§” is the Hankel function of order zero which is 
bounded as ro. It is readily verified that this solution differs from that given by 
(9.9) only by a constant multiplier; one can do so, for example, by using the well- 
known identities involving the derivatives of Bessel functions of different orders. 
Once y is determined we may write (9.8) in the form 


a/9 a er 
— (— - i)e = Ai—— Ho (tkr), A arbitrary. (9.11) 
Oxdy 


This means that our function ¢, if it exists, must satisfy (9.11) as well as (9.2). By 
integration of (9.11) it turns out that we are able to determine ¢ explicitly without 
great difficulty on account of the simple form of the left hand side of (9.11).7° This we 
proceed to do. 

Integration of both sides of (9.11) with respect to x leads to 


a a 
(= - 1)¢ = Ai— Hy (ikr) + g(y), (9.12) 
oy oy 


in which g(y) is an arbitrary function. But g(y) must satisfy (9.2), since all other 
terms in (9.12) satisfy it. Hence d*g/dy?—k’g=0. In addition g(0)=0, since the 
other terms in (9.12) vanish for y=0 because of (9.3) and the fact that 0/dyH\ 
=ik(y/r)dH$/dr. Finally, g(y) is bounded as y>— because of condition (9.6) and 
the fact that 0H! /dy tends to zero as r>~. The function g(y) is therefore readily 
seen to be identically zero. By integration of (9.12) we obtain (after setting g(y) =0): 


‘ 

g=A ie f ote [Ho (ik\/x? + t?) |dt + B(x)ev. (9.13) 
+o 

The function B(x) and the real constant A are arbitrary. The integral converges, 

since 0/dt(H{”) dies out exponentially as t— ». 

We shall see that two solutions ¢;(x, y) and ¢2(x, y) satisfying all conditions of 
our problem can be obtained from (9.13) by taking A =0 in one case and A #0 in the 
other case, and that these solutions will be 90° “out of phase” at ~. (This is exactly 
analogous to the behavior of the solutions in our previous cases.) 

Consider first the case A = 0. The function ¢ given by (9.13) satisfies (9.2) only if 

d? B(x) 
——— + (1 — k*)B(x) = 0. (9.14) 


dx? 





29 It can now be seen how the differential equation corresponding to (9.11) could be obtained for the 
case of waves coming in toward shore over a uniformly sloping beach at inclination 2/2n to the horizontal: 
The left hand side would be a differential operator on ¢ of order 2n, which could be obtained from (5.9) 
by going over to real operators. The right hand side would be essentially the Hankel function Hy) (ikr) 
of order 2n multiplied by sin 2”@. In fact, the entire discussion of sec. 5) for the angles x/2n could be gen- 
eralized to yield three-dimensional progressing wave solutions by proceeding in the manner of the present 


section, 
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It is important to recall that k?<1. The boundary condition g,=0 for x =0 requires 
that B,(0)=0. The condition g,—g=0 for y=0 is automatically satisfied because of 
(9.12) and g(y)=0. Hence B(x) =A, cos~/1—?x, with A, arbitrary, and the solution 
g(x, y) is 


gi(x, y) = Aye” cos y 1— kx. (9.15) 


This leads to solutions ®, in the form of standing waves,* as follows: 


a f cos mf - 
®,(x, y, 2, t) = Ayette” cosv/1 —. k? x. (9.15’) 
: (sin kz 
for k® <1. If R=1, the solution ®, given by (9.15’) is valid, although it could not be 
obtained by our process since B(x) would be identically zero then. 

As we have already stated, we obtain solutions ¢g2(x, y) from (9.13) for A #0 which 
behave for large x like sin+/1—x rather than like cos+/1—2x, and with these two 
types of solutions progressing waves approaching an arbitrary plane wave at © can 
be constructed by superposition. 

We begin by showing that (9.2) is satisfied for all x >0, y<0 by ¢ as given in 
(9.13) with A #0, provided only that B(x) satisfies (9.14). Since x >0, it is permissible 
to differentiate under the integral sign in (9.13), even though ¢ takes on the value zero 
(since the upper limit y is negative). By differentiating we obtain 


f v 0 0? (1) 
V*9 — k®g = Ai: of et — | -+(1—- a) | dt 
\ ie 0tLdx? 


= 2a) 
dHo d Ho 
dy dy 





\ + B(x) +(1— kB. (9.16) 


Since H§” is a solution of (9.2) the operator (02/dx?—k*) occurring under the integral 
sign can be replaced by —0?/dy? and hence the integral can be written in the form 


“ 03 0 2... 
f «| —-+— |p (ikr)dt. 
0 ot? ot 
We introduce the following notation 
a a 
I n(x, vy) = ef —o Ho (ikr)dt, 
= ot* 


and obtain through two integrations by parts the relation 


| a” 1 Pie ™ 7 of F ~~ 
I(x, y) = | + — | + ef e€ - dt, 
+ 


hie ra “1 0 y m—2 ‘is ot m—2 











in which we have made use of the fact that the boundary terms are zero at the lower 
limit + ©, since all-derivatives of H{ (ikr) tend to zero as r++. The integral of 
interest to us is given obviously by J,—J; and this in turn is given by 


30 The standing wave solutions of this type (but not of the type with a singularity) for beaches sloping 
at angles x/2n were obtained by Hanson [3] by a quite different method. 
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2 (1) 1) (1) 
aH dH. v OH 
do Be Bg eae dean et es of e-! ~dt 
ay? ay a al 
- (4) 2. (1) ee) 

y OH» 0 Ho OH> 

+e f e-t - ikea: tie cegaemin bb es 

a ot dy? dy 


by use of the above relations for J,. Hence the quantity in brackets in (9.16) is 
identically zero—in other words the term containing the integral on the right hand 
side of (9.13) is a solution of (9.2). Hence ¢ is a solution of (9.2) in the case A #0 if 
B(x) satisfies (9.14). Since (9.12) holds and g(y) =0 it follows that the free surface con- 
dition (9.3) is satisfied by ¢ in view of the fact that 0H{ (ikr)/dy =0 for y=0. 

We have still to show that a solution B(x) of (9.14) can be chosen so that g,=0 
for x=0, and that ¢ has the desired behavior for large values of r. Actually, these 
two things go hand in hand, just as in former cases. An integration by parts in (9.13) 
yields the following for ¢g: 

7 onl 26) ——_—_— ma... a y 
g= Aier [ e Ho (thy x? + #)dt+ AiHo (thy x? + y?) + B(a)e , (9.17) 
provided that x>0. It should be recalled that the upper limit y of the integral is 
negative; thus the integrand has a singularity for x =0 since t=0 is included in the 
interval of integration and H{(ikr) is singular for r=0. We shall show that 
limz..0~/dx =0 provided that B,(0) = —2A #0. We have, for x>0 and y <0: 


dp a 0... -<T—= _9o 0... —<———J, 
— = Ate e~*— [Ho (iky/ x? + t?) |dt + Ai— | Ho (tkv/ x? + y?) | + B,(x)e. 
Ox “ Ox Ox 

The second term on the right hand side is readily seen to approach zero as x—0 since 
this term can be written as the product of x and a factor which is bounded for y <0. 
For the same reason it is clear that the only contribution furnished by the integral 
in the limit as x0 arises from a neighborhood of t=0 since the factor x may be 
taken outside of the integral sign. We therefore consider the limit 


- 9 as 
lim f e'— [iHy (ik/@+P)dt, ¢€>0. 
z—0 ‘ Ox 


The function iH{ (ikr) has the following development valid near r=0: 


) 2 | 
iHy (ikr) = — — [Jo(ikr) log r + p(r)| | 
T 


in which p(r) represents a convergent power series containing only even powers of r 
(including a zero power), and J) is the regular Bessel function with the following de- 
velopment 
(kr)? 
22 





Jo(ikr) = 1+ oe 


It follows that 
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(1) 2 x x 
[iHo (ikr)| = —- E Jo(ikr) + J¢ (ikr) | log r + xa(0 | 


Ox Tv 
zt 
=-—- | = JoCit + $kx log r + xa(7) | 
x Lr? 
in which g(r) = (1/r)dp/dr is bounded as x—0 since y <0. The contribution of our in- 
tegral in the limit is therefore easily seen to be given by 
. 2 & x ; 2f" 
lim — - é~'- —— dt = lim — — ———- dt 
2-0 Td~ x? + ? 20 rd, x?+7 
By introducing u«=t/x as new integration variable and passing to the limit we may 
write 


2 € x 2 00 du 
lim — - f —— dt = —- f a 
z0 rd. x*?+Ff edo 1+ 


It therefore follows that lim,..0¢/0x =0 provided that 





BAO) = — 2A. (9.18) 
The function B(x) which satisfies this condition and the differential equation (9.14) is 
2A ; — 
B(x) = -—- sin V1 — k? x. (9.19) 
</1 — k? 


Since H}' (kr) dies out exponentially as r— it follows that the solution ¢ given by 
(9.17) with B(x) defined by (9.19) behaves at like e” sin [(1 —k?)*/2%x J. 

A solution ¢» of our problem which is out of phase with ¢; (cf. (9.15)) is therefore 
given by 


e2x(x, y) = Aa] ier f e Hy (ik\/ x? + Pjdt 


. (1) . SS 2e¥ e at 
+ iy (ikyv/ x? + y?) — =H sin V1 — k? x}, (9.20) 
vi- - 


with A; an arbitrary real constant. A standing wave solution ®, is then given by 


cos kz) 


A 9. 20’) 
sin kz) \ 


By taking appropriate values of k progressing waves tending at © to any arbitrary 
plane wave solution for water of infinite depth can be obtained by forming proper 
linear combinations of solutions of the type (9.15’) and (9.20’). For a progressing wave 
traveling toward shore, for example, we might write 

, V1 — k? . 
&(x, y, 2; 2) = A| gi(x, y) cos ke + r ¢o(x, y) sin kz | cost 
(9.21) 


V1 — k? 
—A Ee y) sin kz — =—— ¢o(x, y) cos ks | sin t 
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in which A; and A? in (9.15) and (9.20) are both taken equal to A. The solution (9.21) 


behaves at ~ like Ae” cos (\/1—#x-+kz+1) as one can readily verify by making use 


of the asymptotic behavior of g(x, y) and g2(x, y).#! 


$(%0) 


Z 


a # 2) 





~J 
0 / Z 3 4 


Fic. 15. Standing wave solution for a vertical cliff (with crests at an angle of 30° to shore). 


The special case k=1 has a certain interest. It corresponds to waves which at 
have their crests at right angles to shore. One readily sees from (9.15) and (9.20) that 
as k—1 the progressing wave solution (9.21) tends to 


$(y, 3; 4) = Ae” cos (2 + 2) (9.22) 


that is, the progressing wave solution for this case is independent of x, is free of a 
singularity at the origin, and the curves of constant phase are straight lines at right 
angles to the shore line. 

The progressing wave solution (9.21) has been discussed numerically for k=1/2, 
i.e. for the case in which the wave crests tend at © toa straight line inclined at 30° 
to the shore line. The function g.(x, 0) is plotted in Fig. 15. With the aid of these val- 


*1 We remark once more that the original space and time variables can be reintroduced simply by 
replacing x, y, z by mx, my, mz and t by at (cf. Sec. 3). 
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ues the contours for ®(x, 0, z; 0), i.e. for the surface of the water at time ¢=0, were 
calculated and are given in Fig. 16. The water surface is shown between a pair of 
successive “nodes” of ®, that is, curves for which @=0. These curves go into the z-axis 
(the shore line) under zero angle, as do all other contour lines. This is seen at once 
from their equation (cf. (9.21) with=0). 
V1i-F 
¢1(x, 0) cos kz + rs ¢2(x, 0) sin kz = n = const. (9.23) 


Since @— © as x0 while g, remains bounded, it is clear that sin kz must approach 


Z 
47 





4) / 2 3 + 


Fic. 16. Level lines for a wave approaching a vertical cliff at an angle. 


zero as x—0 on any such curve. That the contours are all tangent to the z-axis at the 
points z=27n, n an integer, is also readily seen. It is interesting to observe that the 
height of the wave crest is lower at some points near to the cliff than it is at © (where 
the value is minus one), so that a saddle point occurs. We observe also that a contour 
for 7=—1 occurs at the right hand edge of Fig. 16. It is believed that the numerical 
calculations are sufficiently accurate to guarantee the existence of such a contour for 
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n= —1; if so, then this contour is likely to be a closed curve, since n—1 at ©. It may 
be that the wave crest is a ridge with a number of saddle points.” 





APPENDIX I. SOLUTION OF THE COMPLEX DIFFERENTIAL EQUATION IN THE 
GENERAL CASE 


In this appendix we give some of the details of the methods used to obtain explicit 
solutions of the differential equation (5.11) for the analytic function f(z) which satis- 
fies the boundary conditions (5.1) and (5.2). The differential equation (5.11) is 


on 
2 


n i 
[I] (eit >*p —1)-f=A—, Areal. (1.1) 
k=1 


The symbol ]] means, as usual, a continued product, and D=d/dz. The boundary 
conditions are 

RMe(iD — 1)f = 0 for zreal, and (I. 2) 

Re(te“**/2"D)f =O on z= re—*z/2n, (I.3) 


in conformity with (5.1) and (5.2) for a beach sloping at the angle 7/2n. 

In order to obtain all phases at © it is sufficient to obtain solutions f;(z) and f2(z) 
for A =0 and for A #0, respectively. 

We begin with the simpler case A =0. In this case the general solution of (I.1) is 
obviously 


fils) = > cxe***, with (1.4) 
ba] 
with 
By = eft (binti/2) | k= jee (I.5) 


The c;, are arbitrary constants. (The quantities 8, are the same as those used in Sec- 
ton, 5, cf. (5.13).) 

We wish to determine the c, so that the boundary conditions (1.2) and (I.3) are 
satisfied. That this really can be done for all points on these lines is far from obvious 
a priori. To satisfy the condition (I.2) we first write 





n 


H(z) = (iD — 1)f1 = > (Be — 1)exe*** = — Do (e*F/” + 1) exe, 
k=1 


k=1 


and observe that in this relation the coefficient c, may be chosen arbitrarily since 
e'**+1=0. We note also that 8B, =8,_x, in which the bar over a number denotes the 
complex conjugate of the number so that e*#*=@e*n»-* if z is real. In order to ensure 
that ReH(z) vanishes for all real z, it is therefore necessary to choose the constants ¢; 


in such a way that 


(2) = — H(z), for z real. 


# It should be pointed out that we are no more able to decide in the present case than we were in the 
two-dimensional cases whether the waves are reflected back to infinity from the shore, and if so to what 
extent. Our numerical solution was obtained on the assumption that no reflection takes place, which is 
probably not well justified for the case of a vertical cliff, but would be for a beach of small slope. 
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To satisfy the above condition, we equate the pair of terms 
(eitFim 4 A )cpete; = — [erin 4 1]c, pee, = (k= 1, 2,---,n— 1) 


Since e#%* and ek are conjugate for z real, we need only require 


intended eitkin _— 1 
(é rkin +- 1)cx = — [eir(m k)/n + 1|¢, iis or Ch = [ef ] = 


[one a ayo 
Therefore, if the c, are chosen so that 
Cy = 4 tan (rk/2n)é,_x, (k = 1,2,---,m — 1) (1.6) 
the boundary condition (1.2) will be satisfied. 


The condition at the bottom requires a similar analysis of the expression 


n 
K(z) = eit( Inf, = >> eft (m1) 1208 ce 20h for 2 = rea~it/2n 
k=1 


or, upon insertion of the special values of z: 


n 
K (re~i=!2”) es >, eit (2h-D/2n¢, ereit (2k—-14Mm) [2m = Lir). 
k=1 


The real part of L(r) should vanish for all r>0. We observe that ere#t@#—1+™/2n js 
conjugate to eré* (2("+1-)—1+"]/2n” Hence our requirement that L(r) be real leads to 








ett (2k—-1) /2nc, —= — eft (2(n+1—k)—l+n)[2ne |) 


or, as one readily sees 


Gs = fay, for & = 1, 2,.°- =, &, (1.7) 


Thus, in order to satisfy the boundary conditions at the free surface as well as 
at the bottom we must impose on the c; the conditions 


Ce = té,-, tan (rk/2n), R=1,2,---,#-1 (1.8) 


Ck = Cn—k+1) Rm i,2,-*-,%. 


We proceed to show that these relations can be satisfied by a set of values c, which 
are uniquely determined within a real multiplying factor. 
From (1.8) we easily obtain the following recurrence relation by taking conjugates 
and eliminating é, 
wk 


Ca—k = 1Cn—k+1 COt cs k= e. ee ,n— 1, 
2n 


which may also be expressed in the form 
T ar kr 
Cn-k = (1)*c, cot — cot —---cot—» for k=1,2,---,n—1. (1.9) 
2n 2n 2n 


For k=n we have the additional relation 
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C1 = Cn. (1.10) 


If we set k=n—1 in (1.9) the cotangents cancel each other and the relation c¢ = (4)"~'c,, 
results; this relation is easily seen to be compatible with (1.10) only if c, is given as 
follows: 


Cn = re tt(n-l)/4 (2.29) 


in which r is any real number. In order to fix the c, in such a way as to satisfy the 
boundary conditions it is therefore only necessary to choose c, in accordance with 
(1.11) and calculate the remaining quantities by using (1.9). The following somewhat 
more convenient form might also be used to calculate the c,: 


Tv 
Ce = reitl(nt) /4-k/2] cot — - - - cot ————» k= 2,3,---,%, 
2n 2n (1.12) 


Ci = Cn 


the first relation resulting through combining the first of (1.8) with (1.9) and noting 
afterwards that the relation holds for k=n since the product of the cotangents has 
the value one in this case and c, for k =n thus has the value given by (1.11). 

We turn next to the case A #0 and, in fact, set A = 1 for the purpose of the present 
investigation. It is clear that the results for any other value of A are obtained simply 
by multiplication by A. A solution f2(z) of the non-homogeneous equations can be 
obtained without difficulty, though the calculations are somewhat laborious. Instead 
of proceeding constructively we prefer to give the solution and then verify that it 
satisfies all conditions, particularly the boundary conditions. The solution f2(z) is* 


n f zBk e~ 
f(z) = 2. afew f oe dt — rie (I. 13) 


k=1 
In this solution the 8; are given by (1.5) and the a, turn out to be multiples of the cx: 
an = cy/(n — 1)!\/n. (I. 14) 


These values of the a, are of course required in order that f2(z) should satisfy the in- 
homogeneous equation. The path of integration for the complex integrals in (I.13) 
has already been given by Fig. 2 of Sec. 4. The path of integration comes from + « 
along the real axis, then goes along a circular arc with center at the origin (leaving 
the origin to the left), and finally along a ray from the origin to the point 2B,. Since z 
lies in the sector for which —7/2n Sarg 2 $0 and the #; are given by (I.5) it follows 
that 28, always lies in the left half-plane. We consider first the condition (I.2). We 


write 


n i n 1 n 
(iD — 1)fe = >> (iB, — 1)an {Ar} +— Dae = M(t) +—- De 

k=1 S k=1 2 k=l 
in which }A,} hasan obvious significance. Since >>." 14x is real, (as one sees from (1.7) 
and (1.14)) the last term is pure imaginary when 2z is real. We must verify that the a; 
are so chosen that the real part of the remaining terms, M(z), vanishes for z real. 


* Compare with (5.15) of Sec. 5. 
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As with the similar problem with the c, above, we must verify that M(2z) = — M(z) 
for z real. In this case it is the terms 


(iB, — 1)ax{Ax} and (iBr~ — 1)an—z{An—x} 


y 


Z-eXriy 








26, 
Fig. 17 


for which the real parts cancel on addition. By evaluating the residue of £e~‘dt/t 
at the origin we find readily that 


f 7Bk et Bk et 
api;Ant = ap ee f —dt — ries = GO, {emf —di+ ries 
\ Cc +00 t C +20 t 


in which the integral on the right hand side of this equation is taken along the path 
C shown in Fig. 17, while the integral on the left is taken over the path C shown in 
Fig. 2. The origin is kept to the right on C instead of to the left. We observe that 
the exponential terms e*® and e*6»-* are complex conjugates for z real, so that —ie~#* 
and +ie*»-* are complex conjugates, since the same is true for 6; and B,_,. In addi- 
tion, the integrals in {A,} and {A,_,} are complex conjugates, as one sees from the 
above relation: if we replace k by n—k in the integral along C it is clear that it be- 
comes the conjugate of the integral along the path conjugate to C, since the in- 
tegrands are the same but the paths are complex conjugates of each other. Conse- 
quently we have only to verify that 


1 
oo (iB n—t = 1)a,_% 





(iB: _ 1)ax = 


But this is the same as the corresponding relation for the c, given above, and hence 
is satisfied by the a, since the a, are the same as the c, except for a real multiplying 


factor. 
To check the condition (1.3) at the bottom we consider the expression 
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n 


n 
eit (m1) 120) f(z) - p eit /2e-ir/2ng.q,{ Ay} +. glein/2g-in/2n } at 
k=1 


1 


nn 
N(sz) + g-leir/2g—in/2n 2? dk 
1 


evaluated for s=re‘*/2", The last term is pure imaginary for z=re~‘*/?" since ) ta, 
is real. It can be shown that ReN(z) =0 for z=re—**/?" by proceeding in the same way 
as above, except that the terms should be paired in a slightly different manner. In 
fact, the terms 


eit /2e-it/2nB a4) Apt and eft /2e—iel2nQ, y ndng1—k} Anti kt 
are negative conjugates in this case, as one can readily verify. 


APPENDIX II: ASYMPTOTIC BEHAVIOR OF fB¥(e-*/t)dt 


In this appendix we prove a number of assertions made in Sec. 5 with regard to 


2Bk et 
f em (fl, (II. 1) 
+00 t 


the behavior of the integrals 


Be = eit(eintt/9, k= 1,2,+-+ ym (IT. 2) 
when z— in the sector S defined by 
, us 
S:;: O2argz2 — 
2n 
The path of integration is given by Fig. 2 of Sec. 4. 
From (11.2) and the definitions of the sector S it follows that 
a/2 S arg 28% S 3n/2. (11.3) 
It might also be noted that arg 28,=32/2 only for k=n, i.e. for By = —1, and 2 real. 
We shall show that the integral (II.1) behaves for z in S and | z| large as follows: 
( e—7Bk / | vs 
ee a fee em GF, < arg 2B; S 7, 
zBk et By Zz 2 / 
f —di~ (II. 4) 
+2 t e *Bk 1 3m 
| at — —1 — + -> | 7< MES’ 
k z 


\ 


in which the dots refer to terms of higher order in 1/z. 
We consider first the case in which 7/2 Sarg 28,7 and begin by integrating 


twice by parts to obtain 
i et e#Bk oe 2Bk “Bk gt 
— dt = — —— + — +2f — dt. (11.5) 
io: 2B (zBx)? m é 


We shall show that 
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| 2Bk e' | | ek 
| — dt| S ¢|—— 
lyn | 2 


for z in the sector S, with c a 
positive number independent of 
z. Clearly, this would suffice to 


show that our integral behaves 


Y 

Q 
at « like —e~**k/zB,.. ; B 

We consider the integral 

f (—*/t®)dt along the closed path 
indicated in Fig. 17. Since the 
integral over the closed path 
vanishes, it is clear that the be- 


havior of our integral as sz « (4) =x 


_— 





can be reduced to the investiga- R. 
tion of the behavior of [(e~*/t*)dt 
over the circular arc PQR as 
ry— ©. The point P corresponds 
to 28; of course. Upon setting Z=Xtk d 


t=r(cos 6+7 sin 6) the integral 





becomes 


es ae dé 
I = if e~reostg—rising . ____. (¥I.6) 
F Tr 


29216 


For | Z| we have, obviously 


} ae | 
\7|s 2) f Ee -@ 7008819 . 
or? 


e rceosé 6 —rcosd 
| T| <2 —— f d@ = 28» with 6 = arg 2px, 
r? 0 r? 


and we may write 


since e77°°88 < e-°o for 0 <9 <6 Sa. We observe that 
| e~t /(28,)? | = r~*e-" cos 5, f= | 2By |, 


and this establishes (11.4) for 7/2 Sarg 28; S72. 
To establish (11.4) for arg 28,>, we write 


26k t 2Bk e~* 
f — dt = 2ri + ni —- dt, 
Cc + 


20 


by evaluating the residue at the origin. The integral over C (see Fig. 17) can be 
treated in the same fashion as above. The only difference is that [(e~*/#*)dt is taken 
over a circular arc such that 6S$6<0 with —r<65< —7/2. The integral J’ for this 
case corresponding to J in (II.6) is the same as that for J except that 6 is replaced 
by —6. The inequalities for | J’| are thus exactly the same as for | Z| since cos @ is an 
even function and 6 lies in the range <6 S7/2. Thus (II.4) is established in general. 











52 J. J. STOKER : [Vol. V, No. 1 


APPENDIX III: TABLE OF E(z) =e* f° (e~“/u)du 
for z=retbk: B, = eit(k/5+1/2)- b=]. --, 7 


The following table of the function E(z) =e fe~“du/u is an extension of the short 
table of this function which has been prepared by the Mathematical Tables Project 
[8]. The path of integration is that of Fig. 2, but taken in the opposite direction. 

The calculations were begun with fourteen six-decimal place values of the function 
E(z) obtained from the above mentioned source, M.T.P. [8]. [We then computed* 
the rest of the numbers by means of power series expansions along each of the seven 
rays in the second quadrant defined by re***;k=1, ---,7; 18737. |] The values are 
believed to be correct to within one unit in the last figure. 


Real part of E(z) 














| } 
| 1 2 3 4 5 6 7 
te % 
| be = ee SS a eet 
1.00 +.2650 +.1731 +.0657  —.0597 —.2063 —.3784 —.5817 
25 | .1901 1013 —.0023 —.1232 —.2644 —.4304 —.6271 
50 | —.1398 0551 = —.0434 —.1580 —.2012 —.4473 —.6315 
75 | .1046 .0242 —.0688 —.1763 —.3006 —.4449 —.6139 | 
| 
2.00 | +.0790 +.0031 —.0844 —.1849 —.2999 -—.4322 —.5849 | 
2S | .0601 —.0117 —.0939 —.1875 —.2936  —.4140 —.5507 | 
50 | 0457 —.0222 —.0993 —.1864 —.2841 — 3933 —.5149 | 
75 | .0345 —.0296 —.1020 —.1831 —.2730 —.3719 —.4795 | 
3.00 | +.0258 —.0348 —.1029 —.1784 —.2612 —.3507  —.4457 
25 | .0190 — .0385 — .1026 — .1731 —.2494 —.3304 — .4141 
.50 0135 —.0411 —.1015 —.1673 —.2377 —.3112 —.3849 
75 | 0090 —.0428 —.0998 —.1615 —.2266 —.2932 —.3582 
} 
4.00 | +.005¢ —.0439 -—.0978 —.1556 —.2159 —.2766 —.3339 | 
a? 0025 —.0445 -—.0956 —.1499 —.2059 —.2612 -—.3119 | 
.50 .0001 —.0448 —.0933 -—.1444 —.1965  —.2470 —.2920 | 
.75 | —.0019 —.0448  —.0909 —.1391 —.1877 —.2340 —.2739 
} 
| 5.00 | —.0036 —.0446 —.0885 —.1340 —.1794 —.2221 — .2576 
1.25 | —.0050 —.0443 —.0861 —.1292 —.1718 —.2111 —.2428 
.50 | — .0061 —.0438 —.0838 —.1246 —.1646 —.2010 —.2295 
75 | —.0071 —.0433 —.0815 —.1203 —.1579 —.1916  —.2173 
6.00 | —.0079 —.0427 —.0793 —.1162 —.1517 —.1830 —.2062 
| .25 | —.0085 —.0421 —.0771 —.1123 —.1459 —.1751  —.1962 
| .50 | —.0091 — .0414 — .0750 — .1087 — .1404 — .1678 — .1870 
| 75 | —.0095 —.0407 —.0731 —.1052 —.1353  —.1610  —.1785 
| 7.00 | —.0099 -—.0400 —-—.0711 —.1019 —.1306 —.1546 —.1708 





* See Sec. 6 for a description of the procedure employed. 
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Imaginary part of E(z) 











| k | 
| 2 3 4 5 6 7 
|r 
; bE eT : —— es 
| 4.00 | —-.609 -—.7965 -—.68100 —.9299 -.9995 <1.0098 1.129 
| .25 ~.6109 —.6723 —.7294 —.7811 —.8263 —.8633 —.8906 
| .50 —.5411 —.5908 —.6345 —.6705 —.6969 —.7110 —.7092 
| 98 —:4849 —.5253 —.5586 —.5826 —.5946 —.5913 —.5675 
2.00 ~.4366 -.4718 -.0688 <-.515 -.9197 -08 -.8 
| 25 — 3009 —.4272 —.4458 —.4532 —.4463 —.4204 —.3691 
| (so | —.3672 —.3897 —.4031 —.4050 —.3919 —.3592 —.3002 
| 75 ~ 3301 —.3577 —.3670 —.3645 —.3469  —.3004  —.2456 
| 3.00 —.3148 —.3301 — .3361 — .3303 —.3094 —.2688 —.2022 
| 25 —-.2935 —.3062 -.305 <—-.3011 —.2778 —.239838 —.166 
| .50 ~.2748  —.2853  —.2864 —.27614 —.2511 —.2075 —.1399 
| .75 2983 ~.28 -.2 -.2363 -—.205 -.08 —.1177 
| 4.00 | —.2435 —.2504 —.2484  —.2354  —.2087 —.1649 —.0998 
> oe —.2302 —.2358 —.2327 —.2188 —.1917 —.1485 —.0853 
| .80 | —.2183  -.2227  —.2187 = —.2042, — 1770 1345 — .0735 
| .75 ~.2074 —.2109 —.2061 —.1912 —.1641 —.1225 —.0639 
| 5.00 | —.1976 —.2003 —.1948 —.1796 —.1528 —.1123 —.0559 
| 9s — 1886 —.1906 —.1846  —.1693 —.1427 —.1034 —.0494 
50 ~ 1803  —.1817  —.1753 —.1599 —.1338 —.0956 —.0440 
15 —~.1728  —.1736 —.1669 —.1515  —.1259 —.0889 .0395 
6.00 | —.1658 —.1661 —.1992 —.14399 —.1188 —.0829 —.0387 
25 — 1593 —.1593 —.1521 —.1369 " —.1124  —.0777. —.0325 
| so | —.1534 —.1529 —.1457 —.1306 —.1066 —.0730 —.0297 
15 —.1478 —.1471 —.1397 —.1247 —.1013 —.0688 —.0274 
700 | —.1426 —.14146 —.1341 —.1194 —.0965  —.0651 —.0254 





TABLE oF e*F 





k Real Imaginary 

1 — .20791 + .97815 
2 — .40674 + .91355 
3 — .58779 + .80902 
4 — .74315 + .66913 

5 — .86603 + .50000 
| 6 — .95106 + .30902 
7 — .99452 + .10453 
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ON BENDING OF ELASTIC PLATES* 


BY 
ERIC REISSNER 
Massachusetts Institute of Technology 


1. Introduction. In two earlier publications':* the author has considered the theory 
of bending of thin elastic plates with reference to the question of the boundary con- 
ditions which may be prescribed along the edges of a plate. The principal result of 
this work was a new system of differential equations for the deformations and stresses 
in thin plates. With this system of equations it is possible and necessary to satisfy 
three boundary conditions along the edges of a plate instead of the two conditions 
which Kirchhoff has first established for the classical theory. 

The physical basis of these results was recognition of the fact that omission of 
the strain energy of the transverse shears is responsible for the contraction of the 
three physical boundary conditions into two conditions,** and that the problem can 
be treated without this omission. 

While the subject is of interest from the point of view of the general theory of 
elasticity,** it is also of some practical importance, in particular with regard to the 
problem of stress concentration at the edge of holes in transversely bent plates. For 
such problems the classical theory leads to results which are not in accordance with 
experiment as soon as the diameter of the hole becomes so small as to be of the order 
of magnitude of the plate thickness,':* while the new equations which take transverse 
shear deformation into account lead to results which are substantially in agreement 
with experiment.’ 

The main purpose of the present paper is to give an account of the author’s earlier 
derivations? in simpler and more general form. While previously an isotropic homo- 
geneous material was assumed, plates of homogeneous or non-homogeneous construc- 
tion are now considered, with elastic properties which in the direction perpendicu- 
lar to the plane of the plate are different from the elastic properties in directions 
parallel to the plane of the plate. 

As a further example of application of the present system of equations, we treat 
the bending of a cantilever plate due to a terminal transverse load. For the homo- 
geneous plate our result represents a minimum energy approximation to St. Venant’s 


* Received Aug. 7, 1946. 

1 E. Reissner, J. Math. Phys. 23, 184-191 (1944), 

2 E. Reissner, J. Appl. Mech. 12, A68-A77 (1945). 
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ing couple and vanishing twisting couple. The two Kirchhoff conditions which take their places are vanish- 
ing bending couple and vanishing of the sum of transverse force and edgewise rate of change of twisting 
couple. 
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solution, while for the non-homogeneous (sandwich) plate the problem appears not 
to have been discussed previously. 

As before, the results are obtained by an application of the basic minimum prin- 
ciple for the stresses and the Lagrangian multiplier method is used to obtain approxi- 
mate stress strain relations. The discussion of the significance of the Lagrange multi- 
pliers is made more precise compared with that given in the earlier work, in accord- 
ance with comments which have been made.® 

2. Statics and strain energy of plates. Let M, and M, be the bending couples 
HT the twisting couple and V, and V, the transverse shear-stress resultants. Let p 
be the surface load per unit of area (Fig. 1). The equilibrium conditions for an ele ent 
dxdy of the plate are then 


OV. OV, OM, OH 0H OM, : 
: . +p=0, : -—— Pin —-4+4. —V,=0. (1) 
Ox Oy Ox oy Ox Oy 


Equations (1) hold regardless of the way in which the stresses are distributed over 
the thickness of the plate. In terms of the stresses, 


h/2 A/2 he} 2 
M. -f zo ,d2, M, -{ zo dz, H = f ZT zy, 
—h/2 —h/2 —h/2 
h/2 A/2 
2 -f r 222, Vy -f Ty2d2. 
—h/2 —h/2 


Equations (1) are three equations for five unknowns. To obtain further equa- 


a ai i 
AS 


(2) 








—-y— 


























Fic. 1. Infinitesimal elements of a plate in interior and at boundary, showing 
orientation of stress resultants and couples. 


tions, use has to be made of the stress strain relations. This is done here through the 
means of the basic minimum principle for the stresses (Castigliano’s theorem of least 
work) according to which the true state of stress is distinguished from all statically 
correct states of stress by the condition that the complementary energy be a minimum.® 


8 J. N. Goodier, J. Appl. Mech. 13, A251—A252 (1946). 

°E. Trefftz in Handbuch der Physik, J. Springer, Berlin, 1927, vol. 6, p. 73 and Z. angew. Math. 
Mech. 15, 101-108 (1935); I. S. Sokolnikoff and R. D. Specht, Mathematical theory of elasticity, McGraw- 
Hill Book Co., Inc., New York, 1946, pp. 284-287. 
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For a material obeying Hooke’s law, and for given surface stresses or displacements, 
the complementary energy is the difference of the strain energy II, and of the work II, 
which the surface stresses do over that portion of the surface where the displacements 
are prescribed. 

Appropriate expressions for II, and I, are 


1 f 1 2 2 
i. = M,+M,— 2»M,M,+ A! H 
Ii k oe, belie i vt 21 + »)H] 


2 1 2 2 
— = PMs + My) + — (Vet vi} dad, (3) 


n 8 


IT, -$ (M,B, + H.B. + V,w)ds. (4) 


The values of the constants D, C,, and C, depend on the properties of the material 
and on the nature of the stress distribution across the thickness of the plate. Ex- 
amples of their calculation for homogeneous and non-homogeneous plates will be 
given later on. 

The functions 8,, 8, and w are the generalized boundary displacements of the 
problem. As II, measures the work of the boundary stresses it follows that 8, must be 
considered as the angle through which the moment M, turns. A corresponding defini- 
tion holds for 8,. For the same reason the quantity # is to be considered as the 
appropriate measure of the transverse deflection of the plate. The precise meaning 
of B,, 8, and #, in terms of weighted averages of the three components of boundary 
displacement U,, U, and W, will be obtained in the following by equating the work of 
the boundary stresses as given by Eq. (4) to the work of the boundary stresses ac- 
cording to the three-dimensional theory and by reducing the expression of the three- 
dimensional theory to Eq. (4) by introducing the assumed variation of the stresses 
over the thickness of the plate. 

3. Variational derivation of the stress strain relations. To make the complemen- 
tary energy II,—II, a minimum subject to the equations of equilibrium (1), these 
equations are multiplied by Lagrangian multipliers \,, A» and X., respectively, and 
integrated over the plate area. The result is added to II, —I], and the variation of the 


resulting expression is made to vanish: 


"(Mz — My, M, — »M, 2(1 + »)H 
| { 5M, + ———_ 5M, + —————_ 6H 
J (a-—»v)D (1 — »*)D (1 — »*)D 
p 1 . eV. eV, 
~ © 6m. +6M,) +— (VAV. + VAV,) + r.o( Sd na a 
e C, Ox oy 
OM. OH 0H OM, 
4 ra ( + — v.) + a (= +—-— v,)} dxdy 
Ox oy Ox oy ) 
‘aw Pe - 
+ p 1B.5M, + BH, + wiV,}ds = 0. (5) 


Applying Eq. (5) to a rectangular plate and eliminating variations of derivatives by 
integration by parts, we find that the boundary values of the Lagrangian multipliers 


must be 
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loa= wv, Mwm=B6. Kk. = By. 


As Eq. (5) also holds for any part of the plate if the boundary displacements referring 
to this part are identified with the displacements occurring in the actual solution of 
the problem, it follows that the Lagrangian multipliers throughout the plate are re- 
lated to the generalized displacements in the interior of the plate through the equa- 


tions 
Aa = w, A = Bz, Ac = By. (6) 


Introducing Eqs. (6) into (5) and integrating by parts, we obtain the variational 


equation 
M,.—vM Op. M,—vM, 0p 
ff ‘|= ds aie Jou. +{ 2 it a "oar, 
J \LQ—»)D C, ax (i-»)D CC, ay 


21+ »)H OB; OB, Vs Ow 
By ¢ Fe av. 
(1— v?)D dy Ox C Ox 


Ve Ow : 
+|—-—-— -8,|a} A dxdy = 0. (7) 
“het oy 


s 


From (7) follow the generalized stress strain relations of the problem: 
OB. op 1+ p 0p OB, 1+p 
M. = D(- = +7 —+—— r), M, = D(’ = +»—+ r), 
Ox dy C. oy Ox Cy 
(8) 
1-—yp OB. OBy Ow Ve Ow cs 
,——— D(+ =). Pp See ey | ae a ce 
2 Oy Ox Ox C. : oy C, 
The conditions along a boundary f,(x, y) =0 are 
8, = B, or M,=M,, 8, = 8. or H,=H,, w= w or V, = V,. (9) 


Equations (9) are the three boundary conditions appropriate to the present theory 
when displacements or stresses are prescribed. They include the case of a free edge 
(M,=H,=V,=0) and the case of a built-in edge (8,=8,=#=0). Appropriate con- 
ditions for more general edge conditions (such as elastic support) may be derived in 
a similar way. 

The five Eqs. (8) together with the three Eqs. (1) represent a complete system of 
equations for the eight functions V,, Vy, M., M,, H, B:, By, w. When C,=C,=« 
they reduce to the customary equations of plate theory. To obtain the appropriate 
(Kirchhoff) form of the boundary conditions in this limiting case one must, however, 
go back to Eq. (3) and therein make C,= ~ before carrying out the remaining analy- 
sis. 
4. Integration of the system of plate equations. It is possible to transform the sys- 
tem of Eqs. (1) and (8) such that integration in terms of harmonic and “wave” func- 
tions is possible. 

The first of the equations in final form is 
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OV, Vy 
—+— =- >». (10a) 
Ox oy 


Two further equations for V,, V, and w are obtained by introducing the first 
three Eqs. (8) into the last two Eqs. (1) and by observing the remaining Eqs. (8) 
and (1). The result is 


— (yD aVvew 1 1\ dp 

V. —- ————__VV,, = — D—— - (1 +»)p(=—-—) = (10b) 
=. ax 2C, Cr/ dx 
(1—y)D_ aV2w ; 1 1\ dp 

— —V?V, = — D— —(1+»)D{- --=)2, (10c) 
2C, dy 2C, C,/ dy 


where V? =0?/0x?+0?/dy?. 

Once Eqs. (10a) to (10c) are solved, the remaining five quantities M,, M,, H, 
8,, By are found from Eggs. (8) by differentiations only. The first three Eqs. (8) may 
be written in the alternate form 


(0? w 0? w DD OF « i+vp v 
UM, = a | + yp )+a-»- —— + d(-— -=)p, (10d) 
Ox? dy tj. ~~ oe “ine 
0? w 07 w D @V, i+ yp v 
M, = - v( —-+ yp- )+a-n= —- + v(- — -=)>, (10e) 
dy? Ox? C. Oy C. Cs 
0? w 1—vD /ovV, OV, : 
H = — (1 — » D— + — - (= 4+ ——-- ), (10f) 
Oxdy 2 C.\ dy Ox 


where 8, and 8, have been taken from the last two Eqs. (8) and use has been made 
of (10a). 

The system (10a) to (10f) is completed by the last two Eqs. (8). 

The solution of equations (10a) to (10c) requires finding a particular integral for 
the load function # and finding sufficiently general solutions of the homogeneous equa- 
tions. The latter is accomplished, as in the paper quoted in Footnote 2, by satisfying 
the homogeneous equation (10a) by means of a stress function x in terms of which 


Ox Ox 





V.=—) =—-— 11 
ay v a (11) 
With 
aw (12) 
; & 
the homogeneous equations (10b) and (10c) become 
a a a _ esa, 
(x — R2V2x) = ——(DV*w), —(x — k*V*x) = — (DV*w). (13) 
dy Ox x oy 


Since Eqs. (13) are Cauchy-Riemann equations, we have 


DV?w — i(x — k?V*x) = @ + iW = f(x + iy). 
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From x —k?V?x = —y follows 
x="n-¥ (15) 
where y; is the general solution of the “wave” equation (with imaginary velocity of 
propagation) 
v1 — RV, = 0. (16) 


Thus, the stress function x is a combination of a harmonic function y and a wave 
function y;. And if the harmonic contribution to x is taken as the imaginary part of 
a complex function f(x+7y) then DV*w is the corresponding real part. From 


DV*w = (17) 


it follows that, when p=0, w itself is a biharmonic function, just as in the theory 
without transverse shear deformation. 

Some applications of these results to the solution of specific problems for isotropic 
homogeneous plates are to be found in an earlier paper.’ 

5. Homogeneous plates. The values of the constants D, C, and C, in the strain 
energy expression depend on the nature of the plate material. Their determination 
will now be carried out under the assumption that the material of the plate is subject 
to the following system of stress strain relations 








dU 1 ; Ve 
i * 2 oe “= (18a) 
oV 1 Ve 
is = Oy - E (ay — 0s) — E, Oz, (18b) 
0U aV 2(1 + v) 
Y= ry sl ea ae (18c) 
ow 1 , 
¢, = re = zE [o, —vilo,+ o,)|, (18d) 
ow 0U 1 
ae tS a (18e) 
ow oV 1 
Yes = —t+— => ae 


dy dz. G, 


Equations (18) stipulate that the plate is isotropic with respect to directions paral- 
lel to the plane of the plate but has elastic properties in the direction normal to the 
plane of the plate which are different. 

The strain energy for a plate of thickness 4 with the stress strain relations (18) 
is given by 


1 h/2 1 
Il, = —fff \— [o, + oy — 2vo,y + 21+ v) rey] 
2 -nj2 LE 


1 2 1 2 2 
+ E. lo, — 2v.0(02 + oy) | + G. [ree + rel dzdxdy. (19) 
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Equation (19) is reduced to (3) by appropriate assumptions regarding the variation 
of the stresses across the thickness of the plate. It is rational to assume that the bend- 
ing stresses vary linearly over the thickness of the plate, while the transverse shear 
stresses vary parabolically: 


4 


M. 3 My 2 Ss «& 
o.=— ay ity — (20) 
h?/6 h/2 h?/6 h/2 








Oy 


~ 2/6 h/2 


Ve zs \? Vy sz \? 
20tefy- eEite  de 
2h/3 h/2 2h/3 h/2 


Equations (20) and (21) satisfy two of the three three-dimensional differential equa- 
tions of equilibrium provided the stress couples and stress resultants satisfy (1). From 
the third of the three-dimensional equilibrium equations, and from the condition that 
the load # is acting on the face z= +h/2, the transverse normal stress g, is obtained: 


= + Z -(- ] (22 
aaa’ ae er =a) ) 


Substituting Eqs. (20) to (22) into Eq. (19) for II,, we find that the integration with 

respect to z may be carried out and that (19) reduces to (3).* The values of the con- 

stants D, C, and C, are found to be 

Eh’ 5 E,h 5 

D=- =» a 
12(1 — »?) 6», 6 

















Gh. (23) 


Equations (23) are introduced into Eqs. (10). There occurs in particular 
1i-»vy D 1 Eh? 
2 Gc M146, 


o(- +y -) 1 Eh? ( + v) ~) 
c. Cj} 10 1—» E, G,)' 


For an isotropic material (E,=E, v,=v, G,=E/2(1+v)) the terms in (24) reduce to 
the values for these quantities which were first obtained in an earlier paper of the 
author? and Eqs. (10a) to (10f) reduce to Eqs. (1) to (V1) of the earlier paper. 

In order to determine the significance of the generalized displacements B,, By and 
w, we write the work of the surface stresses in the form 


h/2 
Il, -$f lonUn + ta2Ue + TnzW |dzds, (25) 
—h/2 


k? = 





(24) 





where U,, U, and W are the actual displacement components of a point of the bound- 
ary. Substituting (20) and (21) into (25), we have 


I oe pg Se Egg ee E (= ) |r} aca (26) 
<f soe h/2 ”  h/6h/2 ° 2h/3 h/2 wan 


Comparison of Eqs. (26) and (4) gives 











* With the exception of a term containing p* which disappears when the variation is carried out and 
which is therefore not evaluated explicitly. 
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_— 6 Pts 6 preg 
Br - U,.— dz, B. nods f U, rage 1z, 
WJ ay." h/2 BS ae hi? 


3 Al? z \2 
v= —{ m1 _ ( - ) [ae 
2h —h/2 h/2 


As Eggs. (26) and (4) hold for any portion of the plate it follows from Eqs. (27) that 
throughout the interior of the plate. 


6? 6 prt os 
B; = — U — dz, By = - f V dz, 
h?J nye = =—h/2 . - nye h/2 . 


3 h/]/2 zZ 2 
w= — w[i _ (=) Jae. 
2hd —_nje h/2 

From Eqs. (28) it is concluded that 8, and By represent quantities which are equiva- 
lent to but not identical with components of change 
of slope of the normal to the undeformed middle 
surface, while w is a weighted average, taken over 
the thickness, of the transverse displacements of 
the points of the plate. Thus, according to the third 
Eq. (28), the present theory leads to approximate val- ; 
ues not for the deflection of the middle surface of the i 
plate but for a weighted average across the thickness h ‘ 
of the deflections of all points of the plate which lie on 
a normal to the middle surface. 

6. Sandwich plates. We consider a composite tye 
plate consisting of a core layer of thickness h and 
of two face layers of thickness ¢. It is assumed that 

j zZ ; 3 | 

t is small compared with h and that the core ma- aor 

oe ; 2 Fic. 2. Infinitesimal element of 
terial is much more flexible than the face material. __, . ene ies 

2 : sandwich plate, showing dimensions 

Under these assumptions the transverse shears are and relevant components of stress. 
predominantly taken by the core plate while the 
bending stresses are primarily taken by the face plates (Fig. 2). 

We take for the strain energy of the composite plate the following expression* 


(27) 





(28) 








pe 











WLLZZZIZZZZ TL 


t 2 
II, = —ff [ozs + Pe — 2vo2,50y,¢ + 2(1 + v) Ty \dxdy 
f 


1 h/2 2 2 
+ —fff [ree,c + Ty2,c\dedxdy, (29) 
2G. ~h/2 


where the subscript f refers to the face layers. The stresses in the face plates are taken 
to be uniform across the thickness and the relations between stresses and couples are 


then, 
* While the assumptions made in what follows should give an accurate picture (within the linear the- 


ory of bending) for combinations such as a foamy core substance and aluminum face plates they will not 
be sufficiently accurate for plates composed for instance of two different kinds of wood. 
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M; M, H 


+ a ¥ Cys = i ? "st? + ; (30) 
h-+ h-+2) (h+2) 








As no stresses @;, @y, Tz, are assumed to be acting through the core material, it 
follows from the differential equations of equilibrium that the transverse shear stresses 
do not vary across the thickness of the core, 

Fs V 
Ue 

Tyz,.¢ =—- (31) 
h h 


Tzz,c ad 


ae 


Substituting Eqs. (30) and (31) into (29), we obtain 
j 1 2 2 : 1 8 2 
i, = [f ——|M,+M,—2»M,M,+ 2(1+»)H?] +——|[V.+ Vy] > dxdy. (32) 
JJ \i(n+i?Ey 2G.h 


Comparison of equation (32) with equation (4) shows that for the sandwich plate 
the constants occurring in the system of differential equations (10) are given in terms 
of the dimensions and elastic properties of the plate as follows 


1 (h+0?E, 


D= —_———- C, = @, C. = hG.; (33) 
2 (1-—»?) 
i1—»vy D t(h+t)?F v i+y (vi(h+ tPE 
k? = — = ios (. - -—*) a oasis (34) 
2 C. 4h(1+>)G . «£ 2h(1 — v)(1 + v)G 


The magnitude of the effect of transverse shear deformation is primarily deter- 
mined by the magnitude of the quantity k. Comparing the first Eq. (24) for the iso- 
tropic homogeneous plate with the first Eq. (34) for the sandwich plate it is seen that 
the effect is of greater importance for the sandwich plate than for the isotropic plate 
whenever 

th E, 21+») 
— om - _ 


2 G, 5 





h?, 


or whenever the ratio E;/G, is greater than the ratio h/t. 

The significance of the Lagrangian multipliers B., 8, and w in the present case is 
determined in the same manner as for the homogeneous plate. One finds here, in- 
stead of Eqs. (28), that in terms of the components of displacement U, V, W, 


AP)--D} 3b@-CD] 


1 h/2 
w= —f Wdz. 
h J _npe 


7. Plate equations in polar coordinates. For the applications to stress concentra- 
tion problems it is convenient to have Eqs. (10) in terms of plane polar coordinates 
rv, 8. Appropriate transformation leads to 
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Orv, OVe 





dailies a = — rp, 36a 
Or 06 P ( ) 
e , 2 dV~ 1 OV? w 1 1\ 0p 
vy, — | Y, - — — — v.| = — D—— = 1+ 90(; -—- ) ) (36b) 
r> 06 r? or as C,,/ or 
: 2 Os 1 
vs — #[ ve -< 4] 
r*> O06 y* 
1 OV?w 1 1\ 1 op 
= — D— — -(1+»v)D ( -— =): —s (36c) 
r 00 3 C,/ r 00 
0" w vy ow vy 0?w OV, 1+ p v 
M,= —D | +— —+—- | + 2k? —- + p(-— —- )p. (36d) 
Or? r or r? 0g? or af Cc. 
1 ow 1 07y 0*w 
M, = - pj. — + — —— + yp — | 
r Or r> 06? or? 
4 | nh ~| 4 o(- +p =) 26 
2k?] — —+— eed Ic 36¢ 
r 06 r Ca CG 4 
: 0/1 dw 1 OV, 0 [Ve 
Hw = — (1-—v)D ( —) + k? {— —+r— (=) (36f) 
Or\r 006 r 06 or\r 
8 Ow Vv, 36 1 dw Vo (36h 
,=-— +2, | ¥=-——+—- (36 
or ' C, (368) - , 0 C, ) 


Equations (36d) to (36f) have been given in the paper quoted in Footnote 2 for 
the case of the isotropic homogeneous plate. Equations (36b) and (36c) have not pre- 
viously been given. They are included in order to facilitate the obtaining of particular 
integrals of the system of equations for load functions of the form p=cos n6f(r). 

Equations (11) which define the stress function x for the solution of the homo- 


geneous equations take on the form 


1 dx Ox 
V»=——) Ve=-—-—-: (37) 
r 06 Or 
Equations (15) and (16) remain unchanged: 
x=vy-y (15) 


where (7, @) is a harmonic function and y; now satisfies the equation 


” #| = gee es =] 0 (16) 
; Or? r Or r2 06? iia 
Also, as before 


DV*w = ¢, (17) 


where now 


o(r, 0) + ip(r, 6) = f(re*). (14) 
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Suitable expressions for ¢, Y and y; are given in Eqs. (42) to (46) of the paper 
quoted in Footnote 2. In the present formulation which includes nonisotropic, non- 
homogeneous plates the quantity h/\/10 in these equations is replaced by the quan- 


tity k defined in Eq. (12) above. 

8. Bending of cantilever plate by terminal transverse load. As an example of the 
application of the formulas of this paper we may treat Saint Venant’s problem of 
flexure of a beam with rectangular cross section.!*:!!:!2 Taking a plate of width 2a 
and length /, held at x =0, acted upon by a force P at x =/ and free of stress along the 
edges y= ta, Saint Venant’s semi-inverse procedure amounts to setting 


M,=V,=,p=0. (38) 
From (10a), 
OV. 
—— = 0, V,=V.A(y), (39) 
Ox 
and from (10c), 
ODV?w , 
og = 0, DV*w = f(x). (40) 
c y 


Introducing (39) and (40) into (10b), we obtain 


d*V, d 
pon ee ee ae | (41) 


dy? dx 


From Eq. (41) it follows, in view of (39) and (40), that 


’ 


y 
V,=C+A cosh ry (42) DV?w = —Cx+t+ B. (43) 


From (43) follows 


.3 2 
Bou = c= + B— + 4(x, 9), (44) 


where (x, y) is a harmonic function which, according to (38) and (10e), is determined 
from the relation 


0° d°*o - 
M,=- EF + (= —Cx+ 8) ] = 0. (45) 
Oy? Ox? 
Evaluation of Eq. (45) leads to the relation 
—1/_z* x* ry ¥ 
Dw =—— (c =f ~) + a = (Ce — B) + Fa + I. (46) 
1-y 6 2 1—v 2 


10S. Timoshenko, Proc. London Math. Soc. (2) 20, 398-407 (1922). 
11S, Timoshenko. Theory of elasticity, McGraw-Hill Book Co., Inc., New York, 1934, pp. 292-298. 
122 A. E. H. Love, loc. cit., pp. 327-346. 
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From Eqs. (42) and (46) follows for the relevant stress couples as defined by equa- 

tions (10e) and (10f), 
M,= (1+ )(Cx — B), (47) H = — Wy+Ak sinh (y/k). (48) 


The five constants of integration A, B, C, F, I are determined by the following five 


conditions 


w= 6,=0 for x=0, y=0 (49) 
M. = 0, | V dy = P for t= '§ (50) 
H=0 for y= +a. (51) 


It is apparent that as in Saint Venant’s theory it is not possible to satisfy the condi- 
tion of complete restraint at the fixed end and also the actual distribution of the termi- 
nal load cannot be prescribed but only its resultant. As a consequence of this the 
solution has general validity only at distances from the ends x =0 and x =/ which are 
at least of the order of magnitude of the width 2a of the plate. 

With 6, from the fourth Eq. (8), Eqs. (49) become 


F C+A | 
if => 0, D = —_ (52) 
C, 


Equations (50) become, with (47) and (42), 
Cl-— B =0, 2aC + 2kA sinh (a/k) = P. (53) 
Equation (51) becomes, with (48), 
— vaC + Ak sinh (a/k) = 0. (54) 


Solving Eqs. (52) to (54) and substituting into Eqs. (42), (46), (47) and (48), we ob- 
tain the following relations for the stresses and deflections 


ss r v (a/k) cosh (y/k) 
V.=— E + (“ —— 1) (55) 
2a 1+ py sinh (a/k) 


PF { 1 (=) 1 = v (2) (: “) 
2aD(i — v*) (2 \1 6 e 2\l l 
k? a/k x) 
+22(14+>—"—_)=h, (56) 


P fsx P —vf[y_ sinh (y/k) - 
M,=— (- — i), (57) H =— — |2 -- —— |}. (58) 
2a\ 1 21+yvLa sinh (a/k) 


< 








g 
| 





Of particular interest is the distribution of shear stress as given by equations (55) 
and (58). For an isotropic plate the results are similar to a known approximate solu- 
tion!®-!! for the beam with rectangular cross section. They reduce in fact to this known 
solution for large values of a/h. 

The maximum transverse shear occurs at the ends y= +a of the plate, 
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ON BENDING OF ELASTIC PLATES 


; ¥ y a/k 
J z.max => — E + TT tate ( caer? “ie 1) |. (59) 
2a 1+ y\ tanh a/k 


The factor in brackets is the correction to the result of elementary beam theory. 
The following table gives numerical values of this factor for the isotropic plate, the 
orthotropic plate and the sandwich plate, in comparison with the exact values for 
the isotropic plate and the known approximate values for the isotropic plate.'®-"! 











TABLE I. Values of Stress Concentration Factor for Transverse Shear in Cantilever Plate. 


























a/k “| .790 | #1.581 | 3.162 | 6.324 | 9.486 | 12.648 
—“_ -f£ .t -rie 2 20 ee" = 
(a/h)x ae 1 1 | 1 | 2 | 3 4 
(a/h) (On | | } 

V 5iG; | 
Eq. (59), =} (|| 1.050 | 1.180 | 1.545 | 2.331 | 3.121 | 3.91 
Eq. (59),=} || 1.040 | 1.144 | 1.436 | 2.065 | 2.682 | 3.33 — 
Appr. »=4 || 1.040 | 1.143 | 1.426 | 1.934 | | 
~— || 4.033. | 1.126 | 1.396 | 1.988 | 2.582 | 3.176 


“yac a 
Exact, v=} 


The magnitude of the shear 7,y parallel to the faces of the plate follows from equa- 
tion (58). rz, is greatest at the points (+4/2, +7) with 7 determined from 


n sinh (a/k) 


cosh — = ————_ 60a) 

os o/h (60a) 
For sufficiently large values of a/k (practically when a/k >3) Eq. (60a) becomes 

n In (a/k) 

—=1-—- (a/k ) (60b) 

a a/k 


and the corresponding value of Hmax is 
v Pp E In (a/k) + “| 


ee aa Ren ai es Ss 6la 
i+yvp 2 


Hex - 


For homogeneous plates Eq. (61a) gives for the shear stress T zy, max =6Hmax/h?, 


2A 4 G v [a a 
——Toanee , y/- gr oe - ma In — a |. (61b) 
3 V/ 10 G. 1+v7Lk k 


The following table contains some values of y and of the factor in brackets in the 
expression for H. 


TABLE II. Location and Magnitude of Maximum Shear Stress Couple and 
Face-Parallel Shear Stress (v=.25). 
































a/k 790 | 1.581 | 3.162 | 6.324 | 12.65 | @ 
n/a | 578 | 504 | .634 ae ae 
sinh (n/k 
Sea 038 29 | .32 | 48 72 1 
a sin (a k) ed aes = Reo ~~ 
Soom — $$ —| == 
(2A /G, | | . 
JS a reves .008 0.2 256 72 2.30 % 
3P G 


exact? G,=G 
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By means of (59) and (58) we may calculate the ratio of maximum shear parallel 
to the plane of the plate and maximum transverse shear. For a homogeneous plate 
we have, in view of (20) and (21), 


ot, 8/2) 4v an sinh n/k a/k es 
ee eee reece oe 
T2:(a@, 0) i+v hla sinh a/k 1 + v\tanh a/k 
and with h4/a from equation (24a) 


T zy(n, h/2) 
T+2(a, 0) 


4 G »p aftn sinh n/k v a/k ie 
(Ae ae ie 
V10V¥ G,1+v klLa sinh a/k 1 + v\tanh a/k 


From equation (62b) follows in particular the limit relation 

















——- tay(n, h/2) 4 G G ; 
lim aaectitine 8" GW a = 1.2664/— (62c) 
a/k—@ T 22(d, 0) Vv 10 G, G, 


which is independent of Poisson’s ratio. Equation (62c) shows the interesting fact 
that, for very thin plates, the horizontal shear may be larger than the transverse shear 
even for isotropic plates. We have confirmed this result for an isotropic plate by an 
exact calculation" in which the factor 1.266 is replaced by a factor 1.342. 

The analogue of Eq. (62a) for sandwich plates is obtained, by means of (30), (31) 


and the first Eq. (34). One finds 


zv.s(”) hG; a sinh n/k a/k - 
: = =— i/ =a ae E = lf + ——( / — — 1) . (63) 
T 22,c(@) i+p 2iG. kLa sinh a/k 1 + »\tanh a/k 


Table II! contains values of the stress ratio as given by equations (62a) and (63) 
for a range of values of a/k and when v=}. 








TABLE III. Values of Ratio of Maximum Horizontal Shear Stress to Maximum Transverse 
Shear Stress for Homogeneous Plates (v=1/4) and for Sandwich Plates. 





- afk” “| 1.581 | 3.162 | 6.324/12.65 | 30 | 100 | »@ 





16/G. max 7; G; max tz | | 
g/ 1G. maxts _ /Gemaxty | gas! 179] .470! .695 | .950| 1.15 


| 


1.266 





iy MAX Te 


ShGy max t- G maxrz, 





Finally, it may be indicated which form the solution assumes in plate theory 
without the transverse shear terms. Equations (55) and (58) become 


P 1 i 
V As -— — —— (64a, b) 


2Zai+typ 2ai+» 
Eq. (57) remains unchanged and Eq. (56) for the deflection loses the terms involv- 
ing k. 
The load P is thus carried in part by transverse shears distributed uniformly across 
the width of the plate and in part by means of concentrated forces at the edges y= +a 
of the plate. As one would expect, no estimate is possible within the frame of the sim- 
pler theory without the transverse shear terms of the actual magnitude of the shear 


stresses which balance the applied load. 








13 E, Reissner and G. B. Thomas, J. Math. Phys. 25, 241-243 (1946). 


SOLUTION OF THE TORSION PROBLEM* 


BY 
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Harvard University 





PUNCH-CARD MACHINE METHODS APPLIED TO THE 





Introduction. Many problems in Engineering and Physics can be reduced to 
boundary value problems, i.e. to problems involving the determination of a function 
which satisfies a given partial differential equation inside a domain and assumes per- 
scribed values on the boundary of this domain. Despite the fact that the solutions of 
many such problems have been known “in principle,” their actual evaluation has re- 
quired such a great amount of computation that it has only been possible to carry out 


calculations of this kind in a few simple cases. 


The creation of modern computational devices such as punch-card machines, 
IBM-Harvard and Bell Telephone Laboratory machines, the ENIAC computer, etc., 
which can carry out rather extensive computations automatically, has changed the 
picture completely. However, these machines exist only to implement theoretical 
methods, and the reduction of these methods to a form whereby the machine can “take 


hold” represents a problem in itself. 


In previous publications the author developed certain theoretical methods (the 
“method of orthogonal functions” and the “method of particular solutions”) for solv- 
ing boundary-value problems. The present paper illustrates the application of or- 
thogonal functions to the solution of Laplace’s equation (0°¢6/dx*) +(0°¢/dy?) =0 


through the use of punch-card machines. 


1. Formulation of a problem in elasticity. The present paper is concerned with a 


method of solving the torsion problem for a bar of uniform cross-section. 


Let x, y, Z denote rectangular coordinates, the axis of Z being perpendicular to the 
cross-section of the beam.' According to Saint Venant the components u, v, w of the 


displacement vector are given by the expressions 
u= —ryZ, v = 7xZ, w = }7G(x, y) 
and the components of the stresses, Xz and Yz by 


Xz = ur[3(0H/dy) — y], Vz = ur|x — 4(0H/dy) |. 


(1.1) 


(1.2) 


Here 7 is the angle of twist per unit length, u is the modulus of rigidity, G(x, y) and 


HH (x, y) are conjugate harmonic functions. 


On the boundary of the cross-section, the function H(x, y) assumes the values 
x?+y?. These conditions determine G and H uniquely within an additive constant 


for G. 


The torsion problem is thus reduced to the “first boundary value problem” of 
potential theory. We shall describe a method of solving it. As an illustrative example, 


we shall determine the function H/ in the case of the domain indicated in Fig. 1. 


* Received June 14, 1946. 


1 This coordinate is denoted by Z rather than z in order that this last symbol can be used to designate 


the complex variable z=x-+1y. 
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If one of the functions, G or H, is known, the determination of the other function 
is very simple. Despite the fact that for the displacement components we need G in 
our example, we determine the function H, since its values are given on the boundary 
and the reader can directly estimate how much our approximate solution differs from 
the required values on the boundary. 

2. The method of orthogonal polynomials. In this paper the method of orthogo- 
nal polynomials will be used for the determination of the function H defined above. 
The mathematical ideas underlying this method have been developed by the author 
in [1],* pp. 57-59 and supplementary Note No. II, [2] and in other papers which are 
listed in these references. However, the necessary computations are very involved; 
the aim of the present paper is to indicate how 
they can be performed with the use of punch- 


| / Kw card machines. 
OA There exist various methods for solving the 


| boundary-value problem of Laplace’s equation. 
The method of orthogonal functions has the fol- 
' % lowing advantages. 

ae iia commer inn wal. 1) One can obtain the solution in a form in 
which the dependence on various parameters is 
evident; in particular one can obtain a solution 
not only for one fixed domain but for a whole 
family of domains which depend upon one of 
several parameters. For instance, in the case of 
the domain indicated in Fig. 1, one’can obtain 
Fic. 1. The cross section B of the bar. formulas involving the radius of curvature at 
the corners and thus investigate how this radius 

















influences the stress distribution. 

2) In order to solve the boundary-value problem, the bulk of the computation, 
i.e. the determination of orthogonal functions, need be performed only once. As soon 
as a set of orthogonal functions for a domain B is known the harmonic function which 
assumes given boundary values or given values of the normal derivative on the bound- 
ary can be easily determined. 

3) The method is very general; it can be easily extended to various other prob- 
lems, e.g. to problems in elasticity as well as to problems involving linear differential 
equations of elliptic type with non-constant coefficients (see Sec. 8). 

Let z=x-+7y, (x, y real), and let ¢‘"+»(z) denote the polynomial 


(**))(z) = Dn+1(2) [E,En+1|-/?, > = 0, 1, Zz ‘tie (2. 1) 


where D,,; and E,4; are the determinants 

| Foo Fo+ ++ Fon—1 1 Foo Fo- ++ Fon—-1 Fon | 
Fae £14°° 325-3 Ss | : Sat bi ; | 
; Eanw=| - y. ieiess ‘ : (2.2a) 


} 
| 
be gs Ros ae ' on 
| Fro Fi. > 2 eae 2” | ive Faja: 7 ae 1 Faual 


Dn+il2) = 


respectively, and 


* Numbers in brackets refer to the bibliography at the end of this paper. 
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Fiyq = ff 2Zdxdy, Z=x-— ty. (2.2b) 
B 


The set of polynomials {¢™(z)}, n=1, 2, 3, - - - , constitutes an orthonormal set 
over the domain B, i.e., 

a 1 if n=m, 
o™goim dxdy = . (2.3) 

B 0 if nAm. 

Now let 
@(")(z) = f o'”(z)dz (2.4) 
0 


and let the real and imaginary parts of 6“ (z) be designated as y, and 6,, respectively, 
so that 
(2) = Yala, y) + 06,(x, 9). (2.5) 

Note that ®‘(z) is a polynomial of degree m in z, and that y, and 6, are polyno- 
mials of degree n in x and y. 

As is well-known, both y, and 6, are harmonic functions. They possess the follow- 
ing property. Let 7(x, y) be a function which is harmonic in B and continuous on the 
boundary C of B. Let s denote the length along the boundary C measured from some 
fixed point, so that the boundary values for H(x, y) may be given in the form 
H(x, y)=f(s) on C. The function H/(x, y) and its conjugate G(x, y) can be expressed 
in terms of prescribed boundary values f(s) and the functions y,, 9, in the following 
form: 


H(x, y) 


I 


at Y| vt, ) f ss)d0.s) — 6,(x, » f sean], (2. 6a) 


G(x, vy) = ¢2+ ie yf f(s)d6,(s) + W(x, yf fovads) | (2.6b) 
y=] Cc Cc 
where c, are constants. The proof of this statement may be found in [1], p. 58. 
In any practical case, the infinite summation (2.6) must be approximated by a 
summation over a sufficiently large number of functions, i.e., (2.6a) must be replaced 


by? 


Hy(x, y) =¢ewt+ D| wee yf f(s)d0,(s) — 0,(x, yf sabe) | (2.7) 


where JN is a positive integer. For increasing N the accuracy of the approximation 
will constantly improve. In the present case, we shall take N=8. 

It is thus seen that the whole problem of solving the boundary-value problem is 
reduced to the performance of the following computations. 

(i) Evaluation of the integrals 

2 In the following some quantities, e.g. 6, and yy, are considered as functions of different variables. 


In passing from one variable to another, new symbols should be introduced, since @, and yy are different 
functions of their respective arguments. For instance 


q) 
6, (s) = 6[x(s), ¥(s)]. 
For the sake of brevity the superscript (or the introduction of new symbols) will be omitted and the func- 
tions will always be denoted by 4,, yy, etc., irrespective of the arguments. 
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ff 23d xdy; (2.8) 
B 


(ii) Computation of the determinants appearing in equation (2.1). 
(iii) Evaluation of the,coefficients 


f f(s)d6,(s), f I(s)dp.(s); (2.9) 


(iv) Determination of the sum appearing in the right hand side of (2.7), i.e. of 
Hy for a sufficiently large N, and the determination of the first partial derivatives 
OHy/dx and 0Hy /dy. 

Using an illustrative example we shall describe in the following how each of these 
operations can be carried out on punch-card machines. 

3. Evaluation of the integrals F,,. In this paper we consider a bar with the cross- 
section shown in Fig. 1. This cross section possesses 8-fold symmetry. 

B is a star-domain; i.e., its boundary C can be represented in polar coordinates 
(p, #) in the form 

p = r(d) (0 <8 S 2rn). £3. 2) 

The values of r(?), 3 =0, 1°, 2°, - - - , 45° are indicated in the second column of 

Table 1. 


TABLE 1: The values of #, r(#), r2(8), Hs(r(d), 8) —cs and of Hs(r(d), 8) on the boundary C of B. 











o r(d) r?(3) H,[r(d), 8] —cs H,[r(a), 3] 
2.722 7.41 —1.309 7.183 
2.423 7.420 —1.295 7.197 
2.728 7.447 —1.249 7.243 
2.737 7.491 —1.174 7.318 
2.749 7.557 —1.070 7.422 
2.764 7.640 — .937 7.555 
2.783 7.745 — .779 7.713 
2.805 7.874 — .595 * 7.897 
2.831 8.020 — .389 8.103 
2.862 8.191 — .163 8.329 
2.896 8.393 -078 8.570 
2.936 8.620 330 8.822 
2.979 8.880 -586 9.078 
3.028 9.169 .840 9.332 
3.083 9.504 1.085 9.577 
3.143 9.878 1.308 9.800 
3.210 10.304 1.498 9.990 
3.276 10.732 1.629 10.121 
3.291 10.837 1.639 10.131 
3.291 10.837 1.610 10.102 
3.273 10.719 1.553 10.045 
3.236 10.476 1.477 9.969 
3.179 10.106 1.388 9.880 
3.122 9.691 1.301 9.793 
3.049 9.303 1.228 9.720 
2.990 8.940 1.163 9.655 
2.941 8.655 1.113 9.605 
2.912 8.486 1.085 9.577 
2. .567 
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The double integral F,,, can be replaced by a single integral. Indeed, 


2 r(g) 
| = f f p™+ ntlei(n—m dd odd 
0 0 
1 


So Te, : r(d) |™+"*2 cos (n — m)ddd 
m+n+ Ay | )| 


+ if [r(d) ]™+"+2 sin (n — moda} . (3.2) 
0 
It follows from the aforementioned symmetry that many of the coefficients Fyn, 
vanish. Indeed, Fn» =0 unless | m —n| =0,4,8, ---.Also, even if F,,, does not vanish, 
it can be shown that the imaginary part vanishes, so that 
1 2 
F,, = F,. = ——— f [r(d) ]™*+"+2 cos (m — n)ddd. (3.3) 
m+n-+2Jo 


If we let m+n+2=s and m—n=q, and take account of the symmetry once again, 
it is easily seen that 


w/4 
Faun = Fam = (8/s) f [r(9) |* cos gddd. (3.4) 
0 


Taking account of the relationship between degrees and radians, and replacing the 
above integral by a summation, we get the approximation: 


Finn = Fam ~ (8/5): (4/180) {ile} + }[r(45°) |* cos (gx/4) 
44 
+ >> [r(v°) }* cos (qre"/180)} ; (3.5) 
v=1 


The above sum can be conveniently evaluated by using the method of “digiting 
without sorting,” see Lorant [4, 5]. The idea underlying this procedure can be ex- 
plained by using the following example. Let us assume that we have to add the follow- 
ing column of products: 

234 X 24 
344 X 32 
342 X 33 
232 X 22. 


In the first column of the multiplicand the number 2 appears in line 1 and line 4. 
Instead of multiplying each number separately, we add the multipliers (24+22 =46) 
and add 46 to itself (i.e., 2*46=46+46). This result is accumulated. Likewise, 
3 appears in the first column of the second and third multiplicands. We there- 
fore add the corresponding multipliers (32+33=65) and then multiply 65 by 3 
(3X65 =65+65+65). This quantity is likewise accumulated. The first column is thus 
completely accounted for; the same procedure is employed for the other columns. 
Thus, the process of multiplication, more time-consuming than addition, has been 
completely replaced by summations, which can be performed on an electric account- 
ing machine. 
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TABLE 2: The values of Fmn. In Table 2-the values of Fun 
: - obtained in our case are given for 
oe ” Fn m=8, n<8. It must be remem- 
sie teak : 0 ns 7 bered that Foo Fan in this case 
, i 122.38 (but not in general). All Frm with 
2 2 737.97 m =8 and n $8 which are not listed 
3 3 5,087.7 are equal to zero. 
: : a 4. Computation of the deter- 
6 6 2,471,500. minants. In order to determine 
7 7 21,041,000. ¢‘(z), defined by Eq. (2.1), we 
8 8 184,060,000. must compute a series of determi- 
4 0 — 169.13 nants whose elements are complex 
5 I 1 —1,523.2 numbers. In this section we indi- 
6 2 —13,843. a 
. 3 ~ 126920. cate a method for evaluating de- 
8 4 —1,174,000. terminants with complex elements 
8 0 —760.2 by the use of punch-card machines. 


Sie hl Tel eae ; We shall again explain the pro- 
cedure by means of an illustrative example, referring for detailed information to 
Lorant [6]. Let the determinant be: 


44, + 1Ayy Gig + tAye G13 + 1A}; | 
G21 + 1A, de2+ iAge dex + 1Ag; | 
431 + 1A31 32 + 1Age 33 + 1A33 
where the quantities a and A are real numbers. 
We shall assume that none of the three complex numbers in the first row is equal 


to zero. We then take the reciprocal of each of the three numbers in the first row. The 
reciprocal of any complex number not equal to zero is given by the formula 


1 a tA ’ 
——_—— = —___ — ——_—__. (4.1) 
a+iA a®?+A* aq?+ A? 
In this way we obtain the three reciprocals 
m+iR, = 1/(ent+idu) (k= 1,2,3) (4.2) 


Now a card is punched for each of the nine elements of the determinant, the entries 
being given in Fig. 2, where m is the number of the row and m the number of the col- 
umn in which the element is located. If a number is negative, we punch a hole in the 
corner above it. See the columns for a, and 7,, where the hole is marked by an X. 
In this case the machines automatically replace addition by subtraction and, in multi- 
plying, punch a hole if one of the factors has a hole. 

Note that the entries in the last two columns will have indices running from 0 to 2 
instead of from 1 to 3, as formerly. Now disregard those cards on which either or 

* The author wishes to thank Dr. Edwin L. Crow for valuable advice and assistance in the perform- 


ance of numerical computations. 
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Fic. 2, The entries on punch cards for the computation of a determinant. 
Every column represents a field on a punch card. 


both of the two subscripts of a* and A* are zero; thus five of the nine cards are dis- 
regarded. From the remaining four cards we set up a new determinant 


* .,* * a 
43, + 1Aqy G12 + 1A, | 





| * oat * i 
| Goi + tAg, doe + tAde 


Repetition of the above procedure reduces this determinant to a single element 


a;**+ iA rn 
The value of the original determinant is then equal to 


(11 + ¢Au1)(a1g + ¢A12)(a13 + 4A13)(0in + 4At1)(Qi2 + GAf2)(Qip + GA). (4.3) 


If we had started with a fourth-order determinant instead of a third-order one, 
the first stage would have reduced it to one of the third order, the second stage to one 


TABLE 3: The values of da, ba, etc. 


























; is bn | te n/n ba/(n—4) | ca/(n—8) 
t || 1.906910-2 = aot | 1.9069 10-3 pa ns 

2 || 9.0407 x 10-2 a ~— 4.5203 X10-2 —_ _ 

3 || 3.6806X10-2 ~ | _ 1.2269 10-2 — _ 

4 || 1.402010 ies is 3.5050 10-3 - sis 

5 || 5.2013x10- | 3.1980x10-2 | = 1.040310- | 3.1980x10-2 on 

6 || 1.8859x10-3 | 2.347210-2 | inn 3.1432X10- | 1.1736X10-% mm 

7 || 6.7240x10-* | 1.261310-2 ~- 9.6057X10-5 | 4.2043x10-3 ‘me 

g || 2.3652x10-~ | 5.900710 | — 2.9566X10-* | 1.4752x10-3 po 

9 || 8 | 1.843710 | 8.1812X10~* | 5.2697x10~+ | 1.8437X10-8 


.2631X10-5 | 2.6348 X10-3 





of the second order, and a third stage would have been necessary to obtain a reduction 
to a single element. From this example it should be quite clear how to proceed in the 


general case. 


Now 


(sz) = ans"! + bys™ + CnZ™-9 + ++ (4.4) 
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and for y, and @, we have the equations 
Vn = (a,/n)r™ cos nd + (b,/n — 4)r"-4 cos (n — 48 + --- (4.5) 


and 
6, = (a,/n)r® sin nd + (b,/n — 4)r™-4 sin (n — 48 +---. (4.6) 


The values of ay, bn, Cn, Gn/n, 6,/(n—4), Cn/(n—8) are given in Table 3. 
5. Determination of the function H(r, 3). In order to determine the harmonic 


90° 














Fic. 3. The exact values of r°(8) along the Fic. 4. The exact values of r*(3) along the 
boundary C of B (first quadrant) and the corre- boundary C of B (first quadrant) and the corre- 
sponding approximate values H,[r(#), 8] indicated sponding approximate values of S[r(8), 8] indi- 
by dashed line. cated by a dotted line. 


function H(r, 3) which assumes the prescribed values, it remains only to determine 
the coefficients. 


A, = f 5 40,s) and B= { ss)a(s) (5.1) 


In our case f(s) =x? +? = [r(8) ]?, where [r(#) ]? is given in column 3 of Table 1. 
The above integral representation for A, can be approximated by the summation 


8r 


A, = — 
180 


{fr(0°) + $[r(45°) ]? + Do [r(u?) ]?[4,(s,) — 0(5,-1)]} (5.2) 


where s, is the point on C having the polar angle of u degrees. 
The expression for B, is exactly the same as (5.2) except that 0, is replaced by y,. 
(It will be noted that instead of summing from 0° to 360° it suffices to sum from 
0° to 45° and multiply by 8; this is not justified in general, but is in this particular case, 
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because of the symmetry of the functions 6, and y, and of the values of r? (=23 +7) 
on C.) 

The summations can be easily performed by using the “method of digiting.” See 
(4, 5]. Because of the symmetry of the specific problem under consideration, most 
of the coefficients vanish. In fact, it can easily be proven, and confirmed by computa- 
tion, that all the A, will equal zero, while only By, Bs, etc., will differ from zero. 
B,and Bs; have been calculated and the values obtained are: By = —4.734, Bg = — 2.343. 

The function H, or more precisely the approximation Hg, obtained by taking only 
the first eight 6, and y,, is now completely determined except for the additive constant 
which we have designated cs. Since the infinite sum has been replaced by a finite 
approximation which we may designate by Hs. it cannot be expected that r? — Hs will 
be absolutely constant. Therefore, the following process is employed in order to ob- 
tain an “average” value of Cg. 

We evaluate Hs on the boundary;i.e., we find >-°_, [A W,(s) —B,6,(s) ] at a number 
of points on the boundary C. We then subtract these values from the correspond- 
ing values of r?, and average the differences. Thus an approximate value of ¢s is ob- 


TABLE 4: The values of H(p, #) in B. 





3 5 1.0 1.5 2.0 2.3 3.0 
0 8.491 8.472 8.389 8.154 7.603 
5 8.491 8.473 8.395 8.177 7.675 
10 8.491 8.477 8.414 8.243 7.874 
15 8.491 8.482 8.442 8.340 8.153 
20 8.492 8.489 8.476 8.453 8.455 
25 8.492 8.496 8.511 8.564 8.727 9.201 
30 8.493 8.502 8.544 8.661 9.936 9.531 
35 8.493 8.507 8.569 8.735 9.074 9.657 
40 8.493 8.511 8.586 8.780 9.147 8.670 
45 8.493 8.512 8 8.795 9 9.661 


.592 


. 169 


tained, and we obtain as the best possible appioximation to H;(r, 3) the quantity 


s 
C3 + y [Aw,(r, 0) — P,6,(r, 3) |. (5.3) 


v=1 


In the example under consideration cs = 8.492. 

In the column 4 of Table 1 the values of Hs [r(d), 3 |—cs and in column 5 those of 
H,[r(2), 3 | are given. The values of Hg(p, 2) in B are given in Table 4. 

6. Methods for improving the approximate solution. The solution obtained repre- 
sents a rather rough approximation.‘ It could naturally be improved by computing 
more orthogonal functions and determining more terms of the development. 

* Since the aim of the present paper is to explain the procedures used rather than to obtain the best 
possible numerical results, the author attempted to avoid long computations and computed only very few 
orthogonal functions. 

In actual applications of the method one will naturally determine more of these functions, in order to 
obtain a much better approximation. 
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In many instances, however, one can save a considerable amount of labor by the 
use of the following simplified procedure: Suppose that there exists a domain, say M, 
which differs only slightly from B, see §1 and Fig. 1., and which includes B+C in 
its interior. Let us assume further that either the function mapping M into the unit 


circle or a complete set of orthogonal functions {f,(z)}, w=1, 2, 3,---, of M are 
known.°® Here we choose for M a circle of radius R=3.295. 
Let (r,, 3,), v=1, 2, - ++, m denote the polar coordinates of m points of the 


boundary C of B, and q, the differences between the desired boundary values and 
those obtained by the use of orthogonal functions, i.e. let 


G = T(r, ,) pe BAfn v,), aa 1, y a ee es (6.1) 


In Table 5, Columns 2, 3 and 6 the values of #,, 7, and q,, respectively are given. 


TABLE 5: The values of r(#), r2(8), Ha(ry, 3»), C(r», 3»), A(%», Jy) and S(r,, Jy) on the boundary C of B. 








dy) | A (ry, dy) | S(ry, dy) 

















v dy ly | ry Hg(ry, dy) | qd c(Ty, 

1 0 2.72 7.41 7.18 | 0.23 | 0.4 —0.02 7.16 
2 10° 2.76 7.64 | 7.56 | 0.10 | 0.5 | 0.07 7.63 
3 20° 2.90 | 8.39 | 8.57 | —0.18 | 1.0 | —0.09 | 8.48 
4 25° 3.00 | 9.00 9.20 | -—0.20 | 2.0 —0.19 9.01 
5 30° 3.14 | 9.88 9.30 | 0.08 | 0.9 | 0.32 | 10.12 
6 35° 3.29 | 10.84 10.13 | —0.S9 | | 0.69 10.82 
7 40° 3.11 9.69 9.79 -0.10 | 0.2 | —0.33 9.46 
8 45° 2.90 8.49 | 9.57 -1.14 | 0.5 —0.98 8.59 
If M is a circle of radius R we introduce the functions 

1 R? — r? 
(R, r3 7,89) = —-——_ —— (6.2) 


2a R? + 7? — 2rR cos (8 — 7) 


where are harmonic for every value of r, and we determine / real constants T,, 


u=1i,2,---+,/,l2m, so that, at the points (r,, 0,), the expression 
I 
A(r, 8) = > T,P(R, 15 1, 9) (6.3) 
p= 


equals q,; i.e. so that 
l 
Dd TyP(R, tui, 3%) =G, v=1,2,--+,m. (6.4) 
p=1 
Remark. The conjugate g(R, 7; r, 3) of p(R, 7; 7, 8) is given by 
rR sin (8 — 7) 
R? — 2rR cos (8 — r) + 7? 





1 
g(R, r37, 8) = — 
T 


5 The function w(z, f) which maps M into the unit circle (taking a point ¢ into the origin) is deter- 
mined, if the set {f,(z)} is known, since according to [1, p. 53] 


wsd= Jf" [Crm ]et[e nol]. 
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The determination of the quantities 7, involves solving a system of m linear equa- 
tions with » variables. This can be performed using punch card machines. The ex- 
pression 

S(r, 3) = H,(r, 8) + A(r, 3) 
can be considered as a second approximation to the desired solution. 

We wish to add that at the points (R, r,), where the difference R—r, becomes 

very small (e.g. for 7, =35° in our example), one has to replace p(R, 7,; 7, 8) by 


(0) ) 7 (1) 
w(R, ‘. é ;r, 0) -f P(R, 7; 7, 3)dr, (6.5) 
7, (0) 
where 7) ri? (7 <r,<7") are suitably chosen quantities. 
We note that the expression (6.3) can be considered as an approximation of the 
integral 
ad 1 sl (R? — r*)h(r)dr 
f P(R, 7; 7, BP) h(1)dr = —{ Sn nna Ee (6.6) 
0 2rJo R*?+ 4r? — 2rR cos (8 — 7) 


where /(r) is some continuous function which at the points 7, assumes the values 
[7,,/(7 —7™) |. In many instances instead of solving Eq. (6.4), we can estimate the 


7 


values of T,, taking T, approximately equal to q,. 





Fic. 5. The values of }°?_,[p(R, 20°+k - 90°; r(8), 8) +p(R, 70°+k -90°; r(v), 3)], 


considered as function of 3, along the boundary C of B. 


In the example under consideration we took for? A(r, &), 


3 
A(r,8) —0.2 >> [p(R, 25° + k-90°; 7, 8) + p(R, 65° + k-90°; 7, 3) | 


k=0 


3 
+ 0.706 >> [p(R, 35° + k-90°; r, 8) + p(R, 55° + k-90°; 7, 3) | 
k=0 


3 
— 1.139 >> p(R, 45° + k-90°; r, 8), (6.7) 


k=0 
R = 3.295. 
7 We note that since r(35°) is almost equal to 3.295, we make an over-simplification by replacing 
0.706|(7s — 76 )w(R, To, ro: 7, 8) + +++ ] by 0.706 [p(R, 35°, 7,8) +--+ J; 
this approximation, although valid approximately for points which are far enough away from 


(R, 35°+k - 90°), (R, 55°+k-90°), k=0, 1, 2, 3, would have to be replaced by more exact expressions in 
the neighborhood of these points. 
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The values of 


c(r, 8) = > [p(R, 25° + k-90°; r, 3) + p(R, 65° + k-90°; 7, 3) | (6.8) 


k=0 


at points (7,,%,) are given in column 7 of Table 5 and drawn in Fig. 5. 

The values of A(r, ®) and S(r, 3) are indicated in columns 8 and 9 of Table 5. 
In the case where M is not a circle, but the mapping function of M into the unit 
circle is known, we can proceed similarly. A very convenient procedure which can be 
applied in this case consists in the use of curvilinear coordinates (in G) which corre- 
spond to polar coordinates in the circle. 

On the other hand, if a complete system of orthogonal functions {f,},4=1, 2, - - 
of the domain M is known, we can construct in M a set of functions P, which possess 
properties similar to the functions p (R, 7,; 7, #). 

The characteristic property of the function p(R, 7; 7,0) is that itis harmonicin M, 
it becomes infinite at the point (r, 3) =(R, 7) and vanishes on the remaining part of 
the boundary of M. 

We denote by P, a system of functions which approximately equal unity on a 
small interval of the boundary and zero on the ren aining part. Using the system }{f,} 
of orthogonal functions, we can detern-ine a set of functions P,, and then proceed as 
before, replacing the p(R, 7,; 7, 8) by the P,’s. 

7. Determination of the stresses. The determination of the stresses requires the 
evaluation of the derivative of H. Our methods yield only functions Hy and S, ap- 
proximating the required function. Inside the domain B the derivatives (0Hy/dx, 
0Hy/dy) and (0S/dx, AS/dy) will in general approximate the corresponding exact 
solution quite satisfactorily. The evaluation of these derivatives may proceed by use 
of punch-card machines, in the manner described previously. Near the boundary or 
on the boundary itself the approximations obtained for the derivatives will, in many 
instances, not be satisfactory, and it is advisable to use some summation method in 
order to obtain more exact values. 

In particular at sharp edges, or points where the radius of curvature of the bound- 
ary is no longer continuous, it is necessary to apply special methods for the evaluation 
of the derivatives. These questions require special treatment, and will be discussed 
in another paper. 

8. Additional remarks. It should be noted that the method described in the 
present paper can be extended without difficulty to the study of thin plates as well 
as to the general problem of elasticity in three dimensions (see [8]). 

Up to the present time, it has been practically impossible to determine Green’s 
function for the equation Au =0 in the case of multiply-connected domains. In a recent 
paper [9] Schiffer has indicated a simple relation between Green's function of La- 
place’s equations and systems of orthonormal functions which permits a compara- 
tively simple determination of Green’s function of such domains. 

Finally by introducing some additional considerations the method of orthogonal 
functions can be extended to a large class of linear partial differential equations in 
two and three variables (see [10], [11]). 
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THE REFLECTION OF AN ELECTROMAGNETIC PLANE WAVE 
BY AN INFINITE SET OF PLATES, II* 


BY 
ALBERT E. HEINS! anp J. F. CARLSON? 
Radiation Laboratory,’ Massachusetts Institute of Technology 


1, Introduction. In Part I* we have calculated rigorously with Fourier transform 
methods, the reflection and transmission coefficients due to the incidence of a plane 
electromagnetic wave on an infinite set of parallel staggered plates. We discussed 
there the case in which there was only one component of the electric field excited; 
that is, the incident electric field was parallel to the edges of the plates. We shall now 
discuss the same geometric structure when it is excited by a plane wave which has 
only a single component of the magnetic field which is parallel to the edges of the 
plates. In this case we shall see that the magnitude of the reflection and transmission 
coefficients are independent of the wavelength and depend only on the angle of stag- 
ger a, and the direction of incidence of the plane wave @. We shall use the Fourier 
transform technique again, and since many of the calculations are parallel to those 
which we did in Part I, we shall only outline the procedure. 

2. Formulation of the problem. We treat here the following problem. A plane 
monochromatic electromagnetic wave whose direction of propagation lies in the plane 

of the paper, is incident upon an infi- 

| nite set of staggered, equally spaced, 

semi-infinite metallic plates of zero 

| thickness and perfect conductivity. 

{/ d These plates extend indefinitely in a 
Py 7 ] ——© direction perpendicular to the plane 
“6 k d of the paper. (See Fig. 1 for a side 

H 1 & i view.) The angle of stagger with re- 
spect to a fixed direction (that of the 

lo cross section of the plates in Fig. 1) 

Fie. 1 is a, while the direction of propaga- 
tion with respect to this fixed direc- 











tion is 8, where a—7<0<a and 6<aS7/2. 

We assume here that the incident wave has only one component of the magnetic 
field, that is, the component which is perpendicular to the plane of the paper. For such 
an excitation, no other components of the magnetic field can be excited. Thus, all 
components of the electric field can be derived from a single component of the mag- 
netic field H.(y, z)=wW(y, z). For this case, the components of the electric field lie in 


* Received April 3, 1946. 

1 Now at the Carnegie Institute of Technology, Pittsburgh, Pa. 

2 Now at Iowa State College, Ames, Iowa. 

3 This paper is based on work done for the Office of Scientific Research and Development under 
Contract OEMsr-262 with the Massachusetts Institute of Technology. 

4 J. F. Carlson and A. E. Heins, The reflection of an electromagnetic plane wave by an infinite set of 
-plates, 1, this Quarterly, 4, 313-329 (1947). Hereafter we shall refer to this paper as I. 
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the plane of the paper and we shall refer to this problem as an “E plane” problem. 
Let us assume as in I, that the time dependence of all field quantities is e~***! 
so that Maxwell’s equations may be written in spatial form as 


VX E = ikH 
VY XH = — tkE, 


where k = 27/X, and ) is the free space wave-length. We then have that ikE, = —dyp/dz 
and ikE,=0y/dy. Upon eliminating E, and E, from the above equations we get the 
two dimensional wave equation 

O*y ay 

—+—+ ky =0 

dy? = az? 


and 


which is to be solved subject to the boundary condition dy/dy=0 on the metal plates, 
since E, is the tangential component of the electric field on these plates. There are 
also conditions at infinity on the function ¥(y, z) which are similar to those which we 
required for $(x, 2) in I. 

We now formulate the equation which expresses the z component of the electric 
field at any point (y, z) in terms of the surface current density on the plates. Following 
a procedure similar to that developed in I, we find 


_ 





+ y’)?]ds’ (2.1) 





¥(y, 2) = ¥is(n 8) — 5 + f I »(2’) 


where y’ =md, g=d cot a and Pine(y, 2 = etk(usin#+z0089) is the incident magnetic field. 
HY [kv/ G —z')?+(y—y’)? | is the free space Green’s function which we have described 
in J. Due to the form in which y and y’ appear, we may write (2.1) as 





v(4 ’ ,z)= a Vine(¥ ’ 2) 


1 
sr >> 7. In(z')Hy [k\/(@ — 2)? + (y — md)*] dz’. 
? mg 


Now 0W/dy, which is proportional to E,(y, 2), is given by the equation 


0 0 
y( vy; = a wkd z) 
Oy Oy 
+7 _f tern VEEP GH mai ee, (2.2) 
y? 
and since E,(y, z) vanishes for y=nd, z>nd cot a,n=0, +1,--- Eg. (2.2) leads to 


an infinite set of inhomogeneous integral equations of the Wiener-Hopf type similar 
to those developed in I. 

Some remarks are now in order about the range of values of d/X which we will 
assume here. In the first place, we suppose that in the parallel plate regions, p(y, 2) 
is proportional to e*** for z sufficiently large and positive. That is, the parallel plate 
region can sustain the so-called principal mode, a mode which cannot occur in the 
H plane case. In order that no higher mode propagate, we must now have that 
0<d/\<4. We further assume that there is only a single reflected wave. This as- 
sumption puts further limitations on d/\ as well as on 6, and we shall discuss them 
when we have obtained the final solution of the problem. 
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3. Fourier transform solution of the integral equation. We have already shown in 
I, Section 3, that the surface current density on any plate can be expressed in terms 
of the surface current density of a given plate, the argument being that we had a 
geometric structure with two types of periodicities. The same argument shows us now 
that 

I n(z a mg) _ To(z)e** ™occst+deing) 

where Jo(z) is the surface current density on the zeroth plate. With this reduction and 
an appropriate change of variables, the integral equation which we have to solve, 


now reads 


oy 





ae : a? a) a 
~~ Z: f To(2’)e**u™ me Hy [kv f{z—2/+(n—m)g}?+(y—md)?]dz’ z>0, (3.1) 
m=—c 0 a 


where now g=d cot a, and y=g cos 6+d sin 6. The y derivatives in (3.1) are under- 
stood to be evaluated at y=nd. Equation (3.1) is an inhomogeneous Wiener-Hopf 


integral equation, which may be rewritten as 


O=7R sin de**#c08° 
_ - ae a Re eee 
re » f To(z')e**#”» — Hy [ky\/(s—2'+ pg)?+(y—nd— pd)? |dz'z>0. (3.2) 
p=—20 0 Oy? 


The dependence of Eq. (3.2) on m is not explicit as the differentiation of 
HS? |kvV/ (z—2/ + bg)? + (y—nd— pd)?| with respect to y will demonstrate. We now 
extend Eq. (3.2) to hold for all z by writing it as 
Pk 0 
F(z) = — Wine(y, 2) 
oy y=0 


where we define 
F(z) =0 z>0O, In(z) =0 2 <0, 


a {= tk sin be'*7 = g > 0, 
Fe Vine( y; 3) : 
me » «36s LO 3 <0. 


C 


I 


We assume as in I, for analytical convenience, that k has a small positive imaginary 


part. 
By arguments which we employed in I, we can show that F(z) is asymptotic to 


eikzcos(2a—-9) for large and negative and therefore it has a Fourier transform 


0 
f e~”2’F(z') dz’ 


which is regular in the upper half of the w plane Jmw > + 9mk cos (2a—6). Since Io(z) 
is asymptotic to e*** for z large and positive, that is, the parallel plate regions are of 
such spacing as to permit the propagation of the principal mode only, the Fourier 
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transform of Jo(z) is regular in the lower half of the w plane {mw < Ymk.> The trans- 
form of OWine(y, 2) /dy| y=o is regular in the lower half of the w plane, ¥mw < Ymk cos 0. 


Finally, the transform of the kernel of Eq. (3.3) is 


Vk? — w? sin d\/k? — w? K_(w) 


— = K(w) = ——— 


cos dy\/ k? — w? — cos (ku — wd) K.(w) 
and it is regular in the strip 
Smk cos (2a — 0) < Jmw < Ymk cos 8. 


Since all of the transforms involved in Eq. (3.3) have a common strip of regularity 
in the w plane, it is thus legitimate to apply the Fourier transform to this equation 
to get 

- k sin 6 Vk? — w? sin d\/k? — w? J(w) (3.4) 
(w) = —_————_ -- sg 3 
w—kcos@ 2[cos dk? — w? — cos (ku — wd) | 

where D(w) and J(w) are the respective transforms of F(z) and Jo(z). 

We are once again led to a problem in factorization of the coefficient of J(w) into 
two functions, one of which is regular in the upper half plane ¥mw> mk cos (2a—8@), 
while the other is regular in the lower half plane ¥mw < Ymk cos 6. The factorization 
of the denominator cos d\/k? —w*—cos(ku—wd) is the same as it was in I, save 
for the fact that p and a have been replaced by uw and d respectively. The numerator 
V/ k? —w? sin dv/ k? —w* may be written in factor form as 


d(k? — w®) T] [V1 = (kd/nn)? + (ited /mn) Je-iwdine 
n=1 


TL [V1 = (kd/n)? — (iwd/nr) Jeivdin, 
n=1 
Clearly (k—w)] [7 [//1—(kd/nx)? +iwd /nx |e-”4/"* has no zeros in the lower half 
plane mw < Qmk, while the factor (k+w)] [2.1 //1—(kd/nx)? — (iwd/nx) Je#4/"= has 
no zeros in the upper half plane Ymw>Qm(—:&). Thus 


Vk? — w? sin d\V/k? — w? 

cos d\/ k? — w? — cos (ku — wd) 
may be factored into two functions, K,(w) and K_(w), where K,(w) is regular in the 
upper half plane Qmw> Qmk cos(2a—6), while K_(w) is regular in the lower half 


plane {im < Ymk cos 6. We have 


K_(2 


& 


)= [ k—w) [] [W1—(kd ‘nn)?+(iwd/nm) jee d/ rex | 
n=1 


i C7 
| (woes Il [A, — iV, Je ea -wot iwod) /2n xt irae) TT [as +i, Jeertieertee—ciet | 
7 -20 


n=1 


5 It has been tacitly assumed that F(s) and Jo(z) are integrable for all finite z. This we can show di- 


rectly from the integral equation. 
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and 
= | k+ Il | 1 —(kd/nx)?—(iwd/ ) Jewwar ( ) | 
— = — Ww ail /n iia 1w / nT) ewd/nxe—x(w 
Kw) d+ ae i i " | 
a_i oo 
| (=o II [An +iW,, Je(tu-wor owe) [San ida) TT [A,— iW, Je(*#-wot iwd)/2nx+ i(hn | 
n=— © n=1 


The A, and YW, are the same functions of w which we encountered in I save for the in- 
terchange of p and a, by uw and d respectively; o,.=k cos 0, ¢.=k cos (2a—8@) while 
x(w)-is an integral function introduced into the decomposition of K(w) to render the 
portions K,(w) and K_(w) algebraic in growth in their respective half planes of regu- 
larity as | w| becomes infinite. 

The decomposition of Eq. (3.4) into two functions, one of which is regular in the 
upper half plane {¥mw>Smk cos (2a—6) and the other of which is regular in the 
lower half plane ¥mw< mk cos 6 follows as before. We now have 

; k sin 0[K,(w) — K,(k cos 6) | 
iets ———— ee. 
w— kcosé 
PEs oR Wire i cle 
2 w— kcosé@ 





Since the left side of (3.4) is regular in the lower half plane {mw < %mk cos @ and the 
right side is regular in the upper half plane {mw> Smk cos (2a—6), while both sides 
are regular in the strip Smk cos (2a—0@) << Xmw<k cos 9, it follows that both sides 
of Eq. (3.5) are equal to an integral function. The final decomposition then gives vs 
k sin 0[K,(w) — K,(k cos 6) ] 


D(w)Ki(w) -— — Integral function, (3.6) 


w— kcos@ 
kK.(kR cos @) sin 0 


— 4](w)K_(w) + Integral function. (3.7) 


w— kcos@ 

We now turn to the evaluation of the integral function of separation and x(w). 
The functional form of x(w) follows from the asymptotic form of K_(w), | w| 0, 
Ymw <0 or of K,(w), | w —«, ¥mw>0. The method of calculating these asymptotic 
forms is the same as in I. We find now that 





K(w)~w?, Ymw <0 | w|—> ©, 


Oe, 





provided x(w) is chosen as —(iwd/r)|(a—x/2 cot a)—In2 sin a]. If we let |w 

Ymw <0, it is clear that the integral function of separation in Eq. (3.7) approaches 

zero in the lower half plane. A similar calculation for Eq. (3.6) shows that the integral 

function approaches zero in the upper half plane | zw | —«, Ymw>0. Since this in- 

tegral function is everywhere bounded and approaches zero for | w| — «, it is identi- 

cally zero. From (3.7) we then get the Fourier transform of I(z) 
2kK+(k cos @) sin 6 


iia ee (3.8) 
K_(w)[w — kos] 


4 


4. Investigation of the far fields. In order to find the amplitudes of the reflected 
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and transmitted waves we must find the asymptotic form of W(y, 2) as | z| —«. This 
can be done by noting that Eq. (2.1) may be written in Fourier integral representa- 


tion as 
1 ee ee 
¥(x, z) = Wine(y, 2) — — dwJ (w)e2t*(ku—wa) n 
4nd c¢ 
cos (y — dn — d)V/k? — w? — e*(*#-“9) cos (dn — y)V/k? — w? 


- (4.1) 


cos Vk? — w?d — cos (ku — wg) 


where the contour C is a path in the strip of regularity §mk cos (2a—6) <<Ymw 
<%mk cos @ and is closed above or below by a semi-circle which passes between the 
poles of the integrand, depending on whether z is greater or less than zero. The in- 
tegration along the semi-circular arc gives no contribution to the integral (4.1) when 
its radius becomes infinite. m is the largest integer in y/d. 

Let us first consider the asymptotic form of ¥(y, z) as z becomes large and posi- 
tive. In the integral in Eq. (4.1), we see that the two poles which give propagating 
mode contributions are at w=o,=k cos 0, and w=k. All other poles give exponential 
factors in z which attenuate for z large and positive. The contour C is, of course, closed 
in the upper half plane. The integral in (4.1) is then 


1 kK.(k cos @) sin de*™#ei(*u—wo)n 
dw il. sli 


(w — k cos 0)\/k? — w? K,(w) sin d\/k? — w? 


2n 
- [cos (y — dn — d)V/k? — w? — e*(*#—“9) cos (dn — y)V/k? — w?] 


= —¢ zcos 6+ysin@) 


ik sin 0K(k cos 0)e***ttkndsing 
EEE Ee (1 oni e*(ku—ka)) 


2kd(1 — cos 6)K,(k) 


+ exponential terms in y and z which are attenuated for z large and positive. 


Thus, for z large and positive, ¥(y, 2) is asymptotic to 
— isin 0K .(k cos O)ettndeiné+sks( 7 i, e*kdsin6) 


oe ee , 


2kdK(k)(1 — cos 6) 





that is, the correct lowest mode functional form in the parallel plate region for “E 
plane” polarization. The magnitude T| of the transmission coefficient is then 
sin a | K,(k cos @) sin [R(u — g)/2] | 


__ sina K is 4.2 
| T| kd(1 — cos 6) | K,(k) | oe 








which reduces to 
sin (a — 6) 


sin (a — 30) cos $0 





|T| = 


provided now & is taken to be purely real. 

We now compute the amplitude of the reflected wave by the same method that 
we employed in I. That is, for y fixed and z large and negative, we close the path by 
a large semi-circular arc in the lower half plane which passes between the poles in 
this half plane. If we allow the radius of this semi-circle to become infinite, the usual 
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arguments show that there is no contribution to the integral from this circular arc. 
The residue calculation then gives us for y fixed, z large and negative, 


ek lysin(2a—6) +z cos(2a@ “OIK(R cos @) sin a] 


¥(y¥, 3) = Vinely, 3) + ————— a Ears = i 
2kK’',[k cos (2a — @)] sin a sin (2a — 0) sin (a — @) 


+ terms which attenuate exponentially for z large and negative. 
With the amplitude of the incident wave taken as unity, the amplitude of the reflected 
wave, that is, the reflection coefficient is 


K(k cos 6) sin 6 








~ 2kK',[k cos (2a — 8)] sin a sin (2a — 8) sin (a — 6) 


which simplifies to 
R = tan 36 cot (a — }6)e* 
when 2 is taken real. ® is the phase angle of the reflection coefficient and is of the 
same functional form as the 0/ —©/ found in I, Section 5, save now for the fact that 
a has been replaced by d. We have finally 

| R| = tan 36 cot (a — 46). 


Conservation of power flow from free space to any parallel plate region gives the fol- 








lowing relation between | R| and | 7]: 
sin (a — 6) 
sage clio bal hol Rol 
sin a 


The values of |R| and |7| clearly satisfy this condition. Finally, we note that the 
condition for a single reflected wave is the same as it was in I, namely 
2d sin a (4.3) 
—_ < ~— e.< 
r cos? [(@ — a)/2] 





Since 0<d/A <3, the inequality (4.3) is not as severe as the comparable relation in I. 
It is illuminating to introduce an angle 7 into the formulas (4.2) and (4.3). This 

angle i is the angle which the direction of propagation of the incident wave makes 

with the normal to the trace of the edges of the parallel plates and is equal to 

ign —a+0. The magnitude of the reflection coefficient then becomes 

sin a — cos i 


| R| = 


aA. 
sin a + cos 7 
while the magnitude of the transmission coefficient becomes 
2 cos 7 


iFji- 


sin a + cos i 


There will be no reflection from this structure if |R| =0 or a=4}n—i, that is, the 
direction of propagation of the incident wave is parallel to surfaces of the parallel 
plates. In this case, of course, | 7| is unity. 





—NOTES— 


A SECOND NOTE ON COMPRESSIBLE FLOW ABOUT 
BODIES OF REVOLUTION* 


By W. R. SEARS (Cornell University) 


In arecent Note! the present author discussed the application of the linear-per- 
turbation theory to the subsonic flow of a compressible fluid past a slender bod of 
revolution. It was shown that the several variants of the theory?:*:+5-* lead to incon- 
sistent results when applied, for example, to the flow about ellipsoids of revolution, 
and it was stated that this ambiguity results from applying the theory to a case that 
is essentially nonlinear. It was implied that all the variants are equally valid so far 
as the first-order theory is concerned. 

In the present Note it will be shown that the various procedures are not, in fact, 
equally acceptable, and that only one of them, in which the boundary condition is 
properly satisfied, can be considered to be correct, even to the first order. 

The boundary-value problem involved is the following: 


Bor: + (1/r)(r¢r)+ = 0, 


Yr dR 
—_—— = — = R(x), say, when r = R(x), (1) 
U+¢; dx 


g—0O as either |x| or r— @, 


Here 6= Ux+¢ is the velocity potential, x and r are the usual cylindrical coordi- 
nates, and 6?=1— M?, a constant defining the Mach number of the undisturbed flow. 
The body of revolution is defined by the function R(x). The perturbation velocities ¢. 
and g, have been neglected in comparison with the free-stream velocity U; similarly, 
yz can be neglected in the first boundary condition in (1). 

It is easily ascertained that solutions I, II, and III of the earlier Note’ do not 
accurately satisfy this boundary condition. These solutions were derived from the 
incompressible-flow solution f(x, 7), which satisfies the following conditions: 

fez + (1/r)(rft)e = 9, 


f(x, S) = U-dS/dx = US'(x), say; (2) 


f—0 aseither |x| —> or ro } 


* Received Sept. 18, 1946. 

1W. R. Sears, On compressible flow about bodies of revolution, Quart. Appl. Math. 4, 191-193 (1946). 

2S. Goldstein and A. D. Young, The linear perturbation theory of compressible flow with applications 
to wind-tunnel interference, British Aero. Res. Com. Reports and Memoranda No. 1909 (1943). 

2H. S. Tsien and L. Lees, The Glauert-Prandtl approximation for subsonic flow of a compressible fluid, 
J. Aero. Sci. 12, 173-187, 202 (1945). 

4H. W. Liepmann and A. Puckett, Introduction to the aerodynamics of compressible fluids, John Wiley 
and Sons, New York, 1946. 

5 R. Sauer, Theoretische Einfiihrung in die Gasdynamik, Springer, Berlin, 1943. Reprinted by Ed- 
wards Bros., Inc., Ann Arbor, 1945. 

* B. Géthert, Ebene und riéumliche Strimung bei hohen Unterschallgeschwindigkeiten, Lilienthal Gesell- 
schaft f. Luftfahrtforschung, Bericht 127, 97-101 (1940). 
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In each solution, the body R(x) is supposed to be related to the different body S(x) 
in such a manner as to satisfy the boundary conditions. The various solutions are 
tabulated below, and in each case is stated the boundary condition that is actually 
satisfied: 

1 
I g(x, r) = f(x, Br), 


¢r(x, S/B) = US'(x); 


II g(x, r) = f(x, Br), 
¢r(x, S/B) = BUS’(x); 
III g(x, r) = f(x/B, 7), 
ox, S) = US"(x/B); 
IV g(x, r) = Af(x, Br), 
v(x, 5/B) = »ABUS"(x); 
V @ = BUx + f(x, Br) = BUx + o(x,7r), say; 


or(x, S/B) = UBS'(x) = (B°U)S'(x)/B. 

It is seen that only Method V, and Method IV with \=1/6?, satisfy in the x, r 
space the correct (linearized) boundary condition for flow about a solid body; i.e., 
that the stream surface has the same slope as the body at the surface of the body. Thus 
the ambiguity mentioned in the Note is eliminated and the formula for the maxi- 
mum superstream velocity ratio on the surface becomes unique: 


Pz _* U 1 - q n 
= = ge fn), (3) 





where is the ratio of maximum diameter to length for the particular family of 
bodies under consideration, and F(m) is the maximum superstream velocity ratio in 
incompressible flow. 

The reason for the confusion mentioned in the earlier Note, and undoubtedly con- 
tributed to by that Note, is not the nonlinear character of F(m), but the fact that ap- 
preciable errors are introduced by use of the approximate boundary conditions shown 
in I-III above, in the three-dimensional case. In the analogous two-dimensional case, 
to be sure, the boundary values of the normal velocity component can be satisfied at 
an approximate location—it is customary to fix them along a straight line, for con- 
venience—without introduction of any error in a first-order theory. Thus the variants, 
which differ only in the location where the slope is fixed, all produce the same result. 
This is not true, in general, in the three-dimensional case. This fact was overlooked 
by various writers*:*.> and by the present author in the earlier Note.! 

The essential difference between the plane case and that of axial symmetry is 
understood by consideration of the distribution of singularities, along the axis, that 
is needed to produce the required flow. In the plane case, both the normal and the 
axial (downstream) velocity components u and 2, are constant, to the first order, 
through a small distance each side of the surface of discontinuity; hence an approxi- 
mation to the position at which the boundary values of v are fixed does not affect the 
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strength of singularities nor the resulting values of u. In the axially-symmetrical case, 
on the other hand, for any source-sink distribution along the axis, ¢, varies as 1/r 
while @, is nearly constant, for small r. Thus an approximation to the radius at which 
the boundary values of @, occur causes appreciable errors in @-.’ 

In conclusion, it now appears, after further study, that the conclusions of the 
earlier Note are in error, and that Method V, or its equivalent, is the only correct 
one for the case of axially-symmetric flow. Method V is that employed by Géthert,® 
and it must be agreed that his objections to the conventional procedure (Method IT), 
which the author previously termed “fancied,” are, in fact, valid.® 


4,.- U 
U 


.20 





°6 ca: «a a wa 7 6 3 1.0 


Fic. 1. The maximum superstream velocity ratio for ellipsoids of revolution in 
compressible flow. M is the free-stream Mach number and 2 is the ratio of maximum 
diameter to length of the ellipsoid. 


For the particular case of ellipsoids of revolution, the results of Method V are 
shown in Figure 1. It should particularly be noted that the correct velocity-ratio 
formula, equation (3), does not permit a universal velocity-correction curve to be 
drawn as function of the Mach number, as the result of Method I would, for example, 
and as in the two-dimensional case. The effect of increasing M is seen to depend inti- 
mately upon m and upon the nature of F(m); hence the curves of Figure 1 cannot be 
expected to apply to bodies other than the ellipsoids of revolution. 


7 In the paper of Goldstein and Young (loc. cit.) this approximation is introduced in equation (21) 
for example, when y is replaced by h/ in integrating to determine y—h. 

8 Note added in proof: In a report recently released [L. Lees, A discussion of the application of the 
Prandtl-Glauert method to subsonic compressible flow over a slender body of revolution, N.A.C.A. Techn. 
Note No. 1127], one of the authors mentioned here has reconsidered the problem and has arrived at con- 
clus‘ons in complete agreement with those of the present note. 
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A NOTE ON STABILITY CALCULATIONS AND TIME LAG* 
By SEYMOUR SHERMAN (University of Chicago) 


In recent technical literature! covering widely different fields, investigations have 
appeared of the zeros of particular exponential sums, one example of which is 


az? + bz + Bze* + ¢. (1) 


The question is: what are the conditions on a, b, 8, and c which are necessary and suf- 
ficient that the real parts of all the roots be negative, thus indicating stability. There 
also have appeared papers in pure mathematics? which discuss similar problems and 
which supply useful techniques for their solutions. It is the purpose of this note to 
indicate how one such technique, which shall be referred to as the Cauchy-Sturm 
method, may be applied to a discussion of the zeros of transcendental expressions 
such as (1). 

Equation (1) arises in the study of control systems with retarded action or time 
lags.* Several attempts have been made to study the zeros of this function and the 
results have not been consistent. Minorsky,‘ in one of his papers, expands the function 


in a power series 


f(z) = az? + bz + Bese“? + 
Bz? z4 B25 Bz8 
+ . ; 


(c+ (6+ B)z + (a — B)2? + _— a oe < 


ll 


He then attempts to approximate the zeros of f(z) by taking zeros of partial sums. 
For nonzero c and 6 we can choose a partial sum of degree n(m23) such that c and 
(—1)"-'8 have opposite signs and so the partial sum 


o3 on 
Zz z 


c+O+ is + fo — Oe + —+--- + (- 
2! (n — 1)! 


* Received March 28, 1946. 

1H, Bateman, The control of an elastic fluid, Bull. Amer. Math. Soc. 51, 601-646 (1945), especially 
pp. 618-626. Further references given there. 

A. Callender, D. R. Hartree, A. Porter, Time-lag in a control system, Trans. Roy. Soc. London (A) 
235, 415-444 (1936). 

A. Callender, and A. B. Stevenson, Proc. Soc. Chem. Industry (Chem. Engg. Group) 18, 108 (1936). 

D. R. Hartree, A. Porter, A. Callender, A. B. Stevenson, Time-lag in a control system, Proc. Roy. Soc. 
London (A) 161, 460-476 (1937). 

N. Minorsky, (i) Control problems, J. Franklin Inst. 232, 519-551 (1941), especially pp. 524-529. 

N. Minorsky, (ii) Self-excited oscillations in dynamical systems possessing retarded actions, J. Appl. 
Mech. 64, A65-A71 (1942). 

F. Reinhardt, Parallelbetrieb von Synchrongeneratoren mit Kraftsmaschinenregelung konstanter Ver- 
zégerungszeit, Wissenschaftl. Veréffentl. Siemens Werke 18, 24—44 (1939). 

2 R. E. Langer, On the zeros of exponential sums and their integrals, Bull. Amer. Math. Soc. 37, 213-239 
(1931). Further bibliography can be found there. 

* For example the roots of (1) determine the stability of the following system with retarded viscous 
term: 

az’’(t) + b2’(t) + B2’(t —1) +c = 0, 

where 2(0) is given and 2’(t)=0, —1<#3s0. 

4N. Minorsky, (i) supra. 
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would have at least one positive zero. This would seem to suggest that dynamical 
systems with retarded viscous terms are necessarily unstable. However, from a theo- 
rem of function theory,’ we know 

For every power series, every point of the circle of convergence is a limit point of zeros 
of partial sums. 
Although function (1) is entire and so has no circle of convergence, the theorem stated 
above shakes our faith in the value of approximating the zeros of function by tak- 
ing the zeros of partial sums of its Taylor series. 

Reinhardt,® in his discussion of the same equation, first considers a particular case 
of equation (1), namely 


227+ 52+ .5ze-* +1 = 0, (2) 


and among the infinite number of zeros of f(z) chooses one which has the largest real 
part and small imaginary part (most unstable and corresponding to low frequency). 
Arguing that this root, corresponding as it does to a low resonant frequency, is physi- 
cally the most significant, Reinhardt studied the zeros of (1) which for other choices 
of parameters could be expanded in a series about the “original” zero. Thus he studied 
one of the infinite number of zeros of (1) intensively and later made approximations 
to the others. He discovered that for some choices of a, b, 8, and c this root had a posi- 
tive real part and for other choices of these parameters this root had a negative real 
part. Thus, the results of Reinhardt and Minorsky are inconsistent. Since both of 
their arguments are approximate a further study is indicated. Minorsky has published 
another analysis’ of this subject which allows for the possibility of stability. 

A method frequently useful for counting the zeros of an analytic function f(z) in 
a simply connected domain D bounded by curve C is Cauchy’s® index theorem: 

If w=f(z) is an analytic function of 2 in a simply connected domain D bounded by 
a closed curve C, f(z) #0, zEC, and z traverses C in a counterclockwise direction, then 
f(z) will traverse a closed curve in the w-plane and the number of zeros of f(z) in D is 
equal to the number of times the w-contour encircles the origin. 

This theorem is at the heart of Nyquist’s® criterion for the stability of amplifiers 
and Routh’s”® stability criterion. An attempt will be made to apply Cauchy’s Theo- 
rem to the zeros of (1). We first note that as z traverses C in a counterclockwise sense 
w may cross the real axis. Let y be the number of times w crossed the real axis in a 
counterclockwise direction relative to the origin (i.e., from quadrant IV to quad- 
rant I or from quadrant II to quadrant III) and let a be the number of times w 
crossed the real axis in a clockwise direction relative to the origin (i.e., from quad- 
rant I to quadrant IV or from quadrant III to quadrant II). The number of zeros of 
f(z) in D is then equal to 1/2(y—a). 

Conditions are sought on a, b, 8, and c (all nonzero) in order that all of the roots of 


w(z) = az* + bz + Bze* +c =0 


5 E. C. Titchmarsh, The Theory of Functions, Oxford University Press, 1939, p. 238. 

6 See (1) supra. 

7N. Minorsky, (ii) supra. 

® H. W. Bode, Network analysis and feedback amplifier design, D. Van Nostrand Co., 1945, chap. 8. 
See Titchmarsh, supra pp. 115-116. 

® See Bode supra. 

10 See E. J. Routh, Dynamics of a system of rigid bodies, Macmillan and Co., 1892, Part II, pp. 191- 
202. 
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have negative real parts. We assume that w(z) has no zeros on the imaginary axis. 
Our domain D is the semicircle 
D: R(z) > 0, |s|<R 
in the z-plane. 
For R(>0) sufficiently large, if | s| =>R and" R(z) 20, then 
| az?| > | Bze~* + bz +c]. 


Since az? €0 for | s| = R, R(z) 20, we have, by arguing from Rouche’s theorem,” that 
w(z) #0, for z in this region. If we choose the R given above, all the zeros of w(z) lying 
in the right half-plane will lie in D and so we apply Cauchy’s Index Theorem to this 
region. The boundary curve C may be broken into two parts 

A; R(z) =0 |s| SR 


and 


B; R(z) >0 | 2] R. 


We consider A for large R. Let z=iy and 
w(iy) = — ay?+c+ Bysin y + i(by + By cos y). 


Note that R(w) is an even function of y and 3(w) is an odd function of y. 

In practical application frequently a>0, c>0. Let us consider the special case 
b=|8| >0. If y=0, then w(0) =c>0. If O<ySR, then 3(w) =y(b+8 cos y) 20 and w 
is in either quadrant I or quadrant II. For large R, w(iR) is in the second quadrant. 
If —Rsy<0, then 3(w)=y(b+8 cos y) $0 and w is in either the third or fourth 
quadrants. Thus, as z traverses A from +2R to —iR, w crosses the real axis once in a 
clockwise direction relative to the origin (from quadrant I to IV). On the other hand, 
for large R, w(z) is dominated by az’ and as z traverses B from —iR to +7R, the net 
number of times that w crosses the real axis is just once in a counterclockwise direc- 
tion relative to the origin (from quadrant IV to I). Since 1/2(y—a) =0 for C, f(z) 
has no zeros in D, therefore all the zeros of f(z) have negative real parts. It should 
be noted that if we remove the restriction of c>0, then" y—a=1-—sgn[c(b+8) |. 

We may rephrase the results as follows: In a one degree of freedom mechanical sys- 
tem with positive mass a, positive spring constant c and positive damping coefficient b, 
and with retarded (unit time lag) coefficient B, if the damping coefficient is greater than 
or equal to the absolute value of the retarded damping coefficient, then the system 1s stable." 

Suppose we relax the restrictions of the preceding two paragraphs but still re- 
quire a>0 and consider (1) on curve C. In order to compute 1/2(y—a) for the line 


11 If gz=x+7y, x, yreal, then 
R(z) =x and Oz) =y. 
2 See Titchmarsh, p. 116. 
13 sen (x) is a real valued function of a real variable defined as follows: 


-i,x 40, 
sgn (x) = 0,x = 0, 
is 8. 


14 This is consistent with Minorsky (ii) p. A69. Note that if 6< —| | <0 (the damping coefficient less 
than or equal to the negative of the absolute value of retarded damping coefficient), then the system is un- 
stable with two zeros of (1) in the right half-plane. 
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segment A and arc B we wish first to note where w crosses the real axis. Because 
R(w(iy)) and 3(w(iy)) are even and odd functions of y respectively, we need only 
consider those positive values of y for which 3(w(ty)) = (6+ 8 cos y) =0. 

If we temporarily ignore the case previously considered where y(b+8 cos y) had 
no positive zeros, but assume merely that 60, then the positive roots of y(b+ 8 cos y) 
will be of the form 


—b 
arc cos (=) + 2mr, m=z=0,1,2,---, 


where we permit arc cos to take two positive values between 0 and 27 (including the 
latter). We number the positive zeros of y(b+ 8 cos y) in order of increasing size 
¥1, Yo, ° * * . Let yw be such a positive zero of b+ 6 cos y that 


(1) 6 sin yw <0 
and 
(2) —ay?+c+ By sin y <0 for!’ y= yy. 


We might, if we so wished, have chosen yy so large that yw <R<ya4:. For the pur- 
poses of our subsequent calculation this would be inconvenient, but since it is con- 
venient for the purpose of the argument, we shall make this assumption during the 
proof. We now calculate y—a for curve C. In other words we have to consider the 
number of times and the direction relative to the origin in which w crosses the real 
axis as 2 traverses C in a counterclockwise direction relative to D. The only values of 
z along A for which w hits the real axis are: 0, +y;;7=0, 1, 2, -- - , M. Let us con- 
sider y;, which is by definition positive. The contribution to y —a due to this crossing 
1S 


d 


— sgn (— ay, +c-+ By; sin y;) sen| }9(b + B cos y) | 


\) 


y= vj 
2 , d 
= — sgn |— ay; +c + By; sin y;]| sgn EF (b + B cos »| 
dy v=; 


= + sgn [(-— ay, +c¢+ By; sin y,;)(6 sin y,) ]. 


Because R(w(iy)) and 3(w(iy)) are, respectively, even and odd functions of y, the 
contribution to y—a, because of the crossing at z= —iy; is also 


sgn [(— ay; +c + By; sin y,)(6 sin y,) J. 
The definition of yy was so arranged that the net contribution to y—a because of the 
crossing of the real axis corresponding to z=0 and those crossings corresponding to 2 
on B totals to 1—sgn[c(b+8) ], as in the case a>0, b=|8| >0, previously consid- 
ered. We now have'® 


If y>(2c/2V/ac—B)>0, then —ay?+c+ vy sin y<0. Thus in many cases the condition 
y >(2c/2\/ac—8)>0 may be used rather than (2). There are other approximations to (2) which might 
prove convenient in different cases. 

16 If 3[w(iy) ] 0 for any real value of y or if 3[w(iy) ]=0 implies that sin y=0, then substitute 0 for 
i. “a sgn } [—ay'+c+by; sin y;][8 sin y;]}. 


ae 
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M 
/ . 2 . . 
(y — a) = 1 — sgn [c(b + 6)] + pe sgn |[— ay; +c + By; sin y;|[B sin y;}}. 
j=1 
Thus a one degree of freedom mechanical system with positive mass a, and nonzero 
spring constant c, damping constant b, retarded (unit time lag) damping constant B, and 
with w(ty) #0, y real, is stable if, and only 2f,"* 


M 
1 — sgn [c(6+ B)| + ‘2. sgn }{|— ay; +c+ By;sin y;|[B sin y;|} =0. (*) 
i=l 


Only minor modifications are required in order to take care of the case a<0 or the 
degenerate cases where one or more of the coefficients is zero. The expression (*) 
may readily be calculated since it involves only four readings from a trigonometric 
table and subsequent evaluations of the sign of quadratic expressions. 

For dynamical systems with a larger number of degrees of freedom or with more 
lag terms, we get higher powers of z or more exponential terms and the application 
is complicated but not hopeless. A machine" of the isograph type should prove helpful 
where extended calculations on complicated systems are being considered. 

Stability calculations are always easier than control calculations. The usual de- 
sign procedure would be to discover a range of a, b, 8, and c corresponding to stable 
responses and then to investigate the detailed response for a few choices of a, }, B, 
and c. Thus we would expect that for any particular choice of a, b, 8, and c the control 
calculation will be more laborious. In the stability calculation of this paper we have 
been, in effect, investigating the relations between the parameters a, }, 8, and c and 
the asymptotic character of the solutions of 

az’’(t) + bz'(t) + Bz’(¢ — 1) + c2(#) = O, 
subject to the boundary conditions: 

2(0) #0, and 2/(4) =0, —1<#80. 
For many design purposes information more specific than the asymptotic character 
of the oscillation might be needed. For instance, after having chosen a, ), B, and ¢ 
so that the transient oscillations are stable (all of the roots of (1) have negative real 


parts), one may be interested in the detailed response 2(t) of the system to a given 
impressed force f(t). This is the control problem. We therefore seek the solution of 


az’’(t) + bz'(t) + B2’(t-— 1) +c = fd), 
subject to the boundary conditions: 
2(0) given, z(t) =0, —1<#s0. 


We consider first the response of the system in the first second. During that time it 
will act like a classical one degree of freedom system with constant mass, viscosity, 
and elastic coefficients a, b, and c, and variable impressed force f(t). The response is 


given by the solution of 


az!"(t) + bz'(t) + ca(t) = f(t), 


17 See T. C. Fry, Some numerical methods for locating roots of polynomials, Quart. Appl. Math. 3, 
89-105 (1945), especially pp. 100-103. 
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where z,(0), given and 2/ (0) =0. The solution of this equation can be found in any 
standard text'* on differential equations, but, depending on the nature of f(t), might 
best be done by numerical integration. In any event, we have 2(t)=2:(t),0<S/31. 
Now during the second second we can again consider our system as a one degree of 
freedom system with the same mass, viscosity, and elastic coefficients, but with an im- 
pressed force which depends on f(t) and velocity at time t—1. Consider the equation 


azz’ (t) + baz (t) + c2(t) = f(t t+ 1) — Bai), 


for OStS1, where 22(0) =2,(1), 2¢ (0) =2/ (1), and 2/ (¢) are derived from the solution 
of the previous equation. Again by standard methods we solve for 22(¢), 0S¢<1. This 
is related to the actual response 2(¢) during the second second by 


2(¢) = 22(¢ — 1), 446 2 


We can continue this process’ for larger values of ¢ if so desired. 


THE CENTER OF SHEAR AND THE CENTER OF TWIST* 
By A. WEINSTEIN (Carnegie Institute of Technology) 


In two recent papers W. R. Osgood! and J. N. Goodier? reconsider the much dis- 
cussed question of the center of shear and center of twist, the former author pointing 
out the disagreement in the literature as to the location of the center of shear. How- 
ever no mention is made of the important paper by P. Cicala* which, together with 
a paper of Trefftz,‘* will form the basis of the following remarks. 

R. V. Southwell’ has clearly pointed out that the two centers, which are intui- 
tively well known to engineers, constitute two different concepts and are not just 
synonyms for the same point. The center of twist is the point at rest in every section 
of a uniform beam subject to a twist by a terminal couple and rigidly clamped at the 
other end. The center of shear (called also flexural center) is the point at which an 
applied shearing force would produce a flexure without torsion. However, neither of 
these points can be explicitly computed, since the displacements of a rigidly clamped 
beam under torsion are not known and, on the other hand, the concept of flexure 
without torsion is still to be exactly defined. Nevertheless Southwell, using Maxwell’s 
reciprocal relations in a summary way, makes plausible the coincidence of both 
centers. 

As Goodier points out, Saint Venant’s theory of torsion and flexure of beams does 


18 E. L. Ince, Ordinary differential equations, Dover Publications, 1944, chap. 6. 

19 Analogous stepwise integration can be found in R. M. Head, Lag determination of altimeter systems, 
J. Aeronaut. Sci. 12, 85-93 (1945) and C. C. Kennedy, Measuring the Coulomb and viscous components 
of friction, Instruments 15, 404-410 (1942). 

* Received May 9, 1946. 

1W. R. Osgood, J. Appl. Mech. 10, A-62—A-64 (1943). 

2 J. N. Goodier, J. Aeron. Sci. 11, 272-280 (1944). 

3 P. Cicala, Atti R. Acc. Sci. Torino 70, 356-371 (1935). 

‘ E, Trefftz, Z. angew. Math. Mech. 15, 220-225 (1935). 

5 R. V. Southwell, An introduction to the theory of elasticity, 2nd ed., 1941, p. 29. 
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not give any definition for either point. Several suggestions have been made concern- 
ing the use of Saint Venant’s theory for this purpose, and it will be shown here that 
a definition can be given, which is based on Saint Venant’s formulae and leads to an 
explicit proof of the coincidence of both centers. 

We shall accept Trefitz’ definition of the center of shear, and refer the reader to 
Trefftz’ paper which also contains a clear exposition of the theory of flexure. It will 
suffice to mention here that Trefftz has been led to his definition by the reciprocity 
law of Maxwell and Betti. 

Let z.be the axis of a uniform beam and x, y the principal directions in the cross 
section S, the origin O being the centroid of the cross section. According to Trefftz, 
the coordinates x» and vr of the center of shear are given by 


1 1 
r= —-— f f yo(x, y)dxdy, Vr = f f xo(x, v)dxdy, (1) 
| 8 I, 8s 


where ¢(x, y) denotes the warping function for the cross section in the theory of torsion 


and 
s Ss 


“ 


are the moments of inertia with respect to the axes of x and y. 

Turning now to the center of twist, the definition which will be given here will 
be based on an idea of Cicala, who however fails to recognize the coincidence of the 
two points as postulated by Southwell. 

According to Saint-Venant, the displacements of the twisted beam are given by 

u=rt(— sy + a+ qz-—Try), 

v=r(sx+6+ rx — pz), (2) 

w = tlo(x, y) +ce+ py — gx, 
where 7 denotes the infinitesimal angle of twist. The last three terms in each formula 
represent a rigid body displacement. We have only six constants a, b, c, p, g, r at our 
disposal, so that a rigid clamping cannot be obtained. However, an approximate 
clamping of the end section z = 0 can be obtained by adopting the following postulates. 

(i) The displacements u and v in the section z=0 are zero: 


u = 0, v=0 for z=0. (3) 


(ii) The mean square of the displacement w in the z direction is a minimum with 
respect to the parameters c, p and g which occur in the last Eq. (2): 


ff wdxdy = minimum. (4) 
Ss 


Denoting this integral by J, we can write in place of (4) 
J(c, p, 7) = minimum with respect to c, p and q. (4’) 


From (2) and (3) we obtain 


y = 0, 





a = 0, 
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so that we have, for any section z=const., 


u=1(q—-y), 0 = 12(x— p). (5) 
We see that, in every section, u and v vanish at the point x =x7, y=yr, where 
xr = p, Vr = q. (6) 


This point will be, by definition, the center of twist. In order to compute p and g, we 
use (4) or (4’), and set the derivatives of J(c, p, g) with respect to c, p and g equal 
to zero. In this way we obtain the following equations. 


ff wdxdy = 0, ff wydxdy = 0, ff wxdxdy = 0, (7) 
¥v¥wv=s5§ S Ss 


Using the last Eqs. (2), namely w=r(y(x, y) +¢+py —qx), and observing that 


f f ydudy = f f adxdy = 0, 
s s 


the origin being the centroid of the cross section, we obtain from (7) the three equa- 
tions 


ff (¢ + c)dxdy = 0, ff (og + py)ydxdy = 0, ff (g — qx)xdxdy = 0. (8) 
Ss Ss S 


The first of these equations yields the value of the constant c, while the second and 
the third gives us the coordinates p=x7 and g=yr of the center of twist, namely 


1 1 
tr = —- ff yo(x, y)dxdy, v= ff xp(x, y)dxdy. 
I. 8 Fy S 


A comparison with (1) shows that the center of shear coincides with the center of 
twist. 

Incidentally, the first Eq. (7), namely [fswdxdy=0, shows that the average dis- 
placement in the direction of the axis of the beam is zero. 

It is interesting to note that both centers, as defined here, are independent of 


Poisson’s ratio. 
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ELASTIC STRESSES DUE TO A SEMI-INFINITE BAND OF 
HYDROSTATIC PRESSURE ACTING OVER A CYLINDRICAL 
HOLE IN AN INFINITE SOLID* 


By O. L. BOWIE (Watertown Arsenal, Watertown, Mass.) 


Recently C. J. Tranter! has presented expressions for the stresses due to a finite 
band of hydrostatic pressure acting over part of the length of a cylindrical hole which 
extends through an infinite elastic solid. By a similar analysis, expressions for the 
stresses due to a semi-infinite band of hydrostatic pressure acting over the cylindrical 
hole can be found. Referring to the analysis by Tranter, it is evident that the only 
change in boundary cenditions is that involving o, at r=a. 

A semi-infinite band of hydrostatic pressure, p, is equivalent to the sum of 


—o <s< @ 
name, f (1) 
2 \ r= 2e 
and 
Perr. —o<s<0 ) 
2 
ty = a, (2) 
ee ceases’ 
2 


The stresses due to the loading (1) are well known,? whereas, the analysis correspond- 
ing to (2) may be carried out in a manner very similar to that used by Tranter. The 
combined stress distributions determined for (1) and (2) yield the following expres- 
sions for the stresses: 


= — £ (=) - “f [apKo(a) Ko(p) + aKo(a)Ki(p) — pKo(o)Ki(a) 


2 \r Tr 


sin (za/a)da 





— {e* + 201 = »)} Kila) Kilo) |- = 


eo 


+2 fi [ako(a)Kilo) — pK o(p)Kx(0)] 
Jo 





Tes cos (za/a)da 


a. D(a) 
p (a pa f” 
w= = (=) = = f [aKo(a)K1(p) + (2v — 1)pKo(p)Ki(a) 
2 a TT J 9 1 
—241- ”)Kx(a) Kilo) | sin (za/a)da 


‘Aa fee ee | 
J [aKo(a)Ko(p) + 2Ko(p) Ki(a) pKi(a) Kile) ]- 





sin (za/a)da 


. 
Il 


where the notation is that used by Tranter. 





* Received Sept. 27, 1946. 

1C, J..Tranter, On the elastic distortion of a cylindrical hole by a localized hydrostatic pressure, Quart. 
Appl. Math. 4, 298 (1946). 

2 A. H. Love, Theory of elasticity, Cambridge Press, 1927, p. 144. 
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It is evident that by a proper superposition of the preceding results, Tranter’s 
expressions for a finite band of pressure may be obtained immediately. A similar 
variation on Rankin’s solution® for the case of a finite band of external pressure acting 
on a solid cylinder will yield expressions for the stresses corresponding to a semi- 
infinite band of external radial pressure acting on a solid cylinder. 


ON A CONFORMAL MAPPING TECHNIQUE* 


By G. F. CARRIER (Brown University) 


1. Introduction. Many physical problems are readily reduced to the problem of 
finding the value of a harmonic function on a smooth closed curve corresponding to 
given boundary conditions and certain interior singularities. For example, it has been 
shown |1] that the problem of finding the flow of an incompressible non-viscous fluid 
past a periodic array of airfoils' may be replaced by that of determining the flow 
within a smooth closed stream-line C with a source-vortex and a sink-vortex at speci- 
fied interior points.? Analogous problems can obviously occur in other physical situa- 
tions which lead to the Laplace equation. 

Here we shall develop a method of mapping a smooth closed curve C conformally 
onto the unit circle so as to carry two arbitrarily specified interior points into the 
points +a. The mapping is continuous within and on C. 

Since the solution requires the integration of a single non-homogeneous Fredholm 
integral equation, we also present a numerical procedure which is useful in solving 
such equations. The over-all procedure seems preferable, in general, to Theodorsen’s 
method of mapping “nearly circular regions” [2], [3], since the equation used here 
has a non-singular kernel in contrast to those in his two simultaneous integral equa- 
tions. 

2. The integral equation. Let us consider the problem of finding that complex po- 
tential F(¢)=¢+iy on the closed curve C# (see Fig. 1) which corresponds to two 
equal and opposite interior logarithmic singularities at the points P and Q and let ¥ 
vanish on C. We may write then, for points on and within C, 


F(¢) = In [(¢ — P)/(E —- OI] +I, (1) 


where f is analytic within C and continuous on C. In principle, no difficulty, would be 
encountered if the unit coefficients of these singularities were replaced by complex 
numbers, but the present form is sufficiently general for our purpose. 

Consider the integral 

3 A. W. Rankin, Shrink-fit stresses and deformations, J. Appl. Mech. 11, A-77 (1944). 

* Received Sept. 13, 1946. 

1 The airfoil shape being arbitrary. 

2 C, in this case, is a non-self-intersecting convex curve with a continuously turning tangent. By 
“smooth” we shall henceforth imply such a curve. 
3 C is any closed, not self-intersecting, curve. 
















102 NOTES [Vol. V, No. 1 


F(¢) 
T (fo) -$ —— df = 0, (2) 
ee fot 


where S is that closed contour (Fig. 1) which contains no singularities. As a conse- 














Fic. 1. A is an arbitrarily chosen origin of the coordinates s, t; the point Z 
may also be taken arbitrarily. 


quence of the fact that the logarithmic terms differ by 277 on opposing sides of the 
cuts, the integral of Eq. (2) may be written in the form 
es ae. a a = 
— wif (So) + J, F(g)diln (¢ — £0); + 207 In ——— = 0, (3) 
$0 Q — $0 
provided the limiting process is performed wherein the radii of the circles about P, Q, 
and ¢ 9, and the width of the cuts, are made to tend to zero. Since Y=0 on C, the 
imaginary parts of the foregoing equation lead to the integral equation 


1 
¢(s) = —f o(t)da + 2 Inr,/r,. (4) 
7 Cc 
The symbols not used previously are defined in Fig. 1. The angle @ is obviously a func- 


tion of s and #, and its increment in Eq. (4) corresponds to an increment dt. This equa- 
tion has been obtained previously by a somewhat less straightforward method [1 | 
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wherein the s.ngularities were taken with complex coefficients.‘ In this case, Eq. (4) 
has terms —T’,8,/m and —I',6,/r. 

A method of solving this integral equation will be presented in Sec. (4). Note 
that for a circle, the kernel is constant (da/0@=4, where =e), and the solution is 
already obtained. A closely analogous situation was previously obtained by Prager 
[4]. 

3. The mapping problem. We now turn to the problem of mapping C on the unit 
circle in the z plane so that P, Q go into the points +a. The mapping is to be regular 
within and on C. The magnitude of the real quantity @ cannot be arbitrarily chosen 
but will follow from the analysis. We first show that such a mapping exists. Two facts 
are well known: (1) there exists a mapping such that the closed curve C goes confor- 
mally into the unit circle in the z; plane and such that P, Q go into the (unspecified) 
points P’, Q’; (2) there exists a bilinear transformation such that the unit circle goes 
into the unit circle and any two interior points P’, Q’, gointo +a (a being determined 
by the values of P’, Q’). Thus the mapping exists. 

Now consider that the integral equation of the foregoing section has been solved 


for $(s). Denote its minimum value by @min. 


Note that for the points z=e*, 
1 
—(z—a)izs— 
a 


F(z) = $2(0) + iO = In -- + K(a) (5) 


1 
(+ a)(s+ :) 
a 


is the solution of the potential problem wherein Y=0 on the unit circle and unit 
singularities are placed at +a. The real number K(a) may be chosen so that the small- 
est value taken by $2(@) is @min.® 

It is evident that when the curve C is mapped into the z plane in the specified 
manner, the potentials ¢, ¢2, in the two planes must have the same value at corre- 
sponding points, i.e. at points 2(¢;) and ¢, Thus, @ must be chosen such that the 
greatest value of @ is equal to the greatest value of ¢2. Once this has been done, the 
mapping function z(¢) along C is implied by the relation 


$2(8) = $(s), (6) 


One may now quickly find 2(¢) along C, and, if it is required, find 2(¢) at any interior 


point by using the Cauchy formula 


The physical problem may now be treated by the conventional methods used 
when conformal mapping is applied. In particular, for the fluid mechanics problem 
mentioned previously, complex singularities may be placed at +a, the imaginary 
coefficient of the singularity corresponding to the point infinitely far downstream in 


‘In [1], the vertex terms have the incorrect sign. 
5 Actually K(a), a, are determined simultaneously by this condition and that appearing in the next 


paragraph. 
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the physical plane being determined by the Joukowski condition.* The velocity dis- 
tribution along the vane profile is determined using only the analytic’ transformations 
which led to C and the function d6/ds. 

If it is only necessary to place one singular point P within the unit circle, it may 
be mapped into the origin. One uses doublets at P and at the origin. The strength of 
the doublet must be adjusted as was a in the preceding problem. 

4. The solution of the integral equation. An elementary numerical procedure 
which resembles the relaxation process conventionally applied to differential equa- 
tions [5] is useful in solving Eq. (4). We choose several points s, (or f,) roughly uni- 
formly spaced along C, compute a(s;, tn) =akn for each k and n, and determine the 
quantities Ajn = (Qk ,n41—Qk,n-1)/27. We now write Eq. (4) in the form 


n rq(Sk) 
and guess the values of $(¢,). Using any convenient arrangement, we compute ¢(s;) 
and compare the computed and guessed results. We then make a better guess for the 
quantities $(¢,) [according to the form of the A,, and to any experience gained in 
previous attempts, computing the change in ¢(s,) | in such a manner that the two sets 
of ¢ values become essentially alike. The accuracy which can be obtained is a function 
of the fineness of the spacing of points in the kernel. 


The usual iteration procedures which are applied to Fredholm equations, or to the 
systems of algebraic equations to which they may be reduced, differ from the foregoing 
in that no freedom of choice is allowed beyond the first approximation. In this 
method, the higher approximations may lead to much more rapid convergence than 
such iterations just as in the relaxation method in the numerical solution of differ- 


ential equations. 
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ON HYPERSONIC SIMILITUDE* 
By WALLACE D. HAYES (North American Aviation, Inc.) 


A recent paper by H. S. Tsien! presents a law of similitude for two-dimensional 
potential hypersonic flows over a slender body. Hypersonic flow is fluid flow for 
which the Mach number is much greater than one. If the transformation 


x= bi, (fa) y= én, (1b) $ = adf(é,n) (1c) 


is made, where 6 and 6 are the length and thickness of the body and a» and M are 
the velocity of sound and the Mach number in the undisturbed stream, the two- 
dimensional potential equation is transformed into 


[1 — (y — 1)Kfe — Bly + Dfalfon = Kee + 2K afew (2) 

where 
K = Mi/b (3) 
is a fundamental similarity parameter. The boundary conditions may be expressed 
fe=f,~=9O at E=—--@, (4a) 
f, = KH’(é) at = H(&), (4b) 


where y=/H-6 defines the body shape. Since the equation is not linearized it is not 
permissible to satisfy the shape boundary condition at 7 =0. Two potential flows with, 
similar bodies and the same value of K are thus given by the same mathematical 
solutions and are similar. The drag and lift coefficients for bodies of similar shape 


may be expressed 


1 : 
Cp = 7) A(K) (5a) 
C A(K) (5b) 
a ——— J pe) 
MM? 


Professor Tsien also demonstrated the analogous law for axially symmetric flow. 

It is the purpose of this note to point out that this similarity law is both much 
simpler in concept and much more general in scope than has been previously indi- 
cated; it is; in fact, applicable to three-dimensional flow with shock waves and 


rotation. 

The two-dimensional hypersonic potential equation expressed by Eq. (2) is 
identical with the one-dimensional non-stationary potential equation in y and ¢ under 
the transformation 


t= Tt (6) 


with (1b) and (1c), provided the replacement 


T = b/Ma (7) 


* Received Jan. 16, 1947. 
1H. S. Tsien, Similarity laws of hypersonic flows, Jour. Math. Phys. 25, 247-251 (1946). 
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is made in the definition of K. Also, the boundary conditions (4a), (4b) are the same 
if now y=H-6 defines the position of the single boundary as a function of time. This 
shows that with the slender body hypersonic assumptions, M> 1, there is no intrinsic 
dependence upon the axial variable from the point of view of an observer stationary 
with respect to the undisturbed fluid and the problem becomes a non-stationary one 
in one fewer space variables. A change of 6 with K kept constant is merely a scale 
transformation in the non-stationary system. 

The difference between this point of view and the actual case is that the disturb- 
ances at two points on the same streamline on the body are assumed to be in phase. 
This difference is negligible if the ratio of the signal time between the two points to 
the time phase difference between the two points is large. Using the fact that in the 
hypersonic flow over a slender body there are appreciable changes in the velocity of 
sound but not in the flow velocity, the ratio of these times is equal to the local Mach 
number, and this parameter being large ensures the validity of the point of view. 
It is clear that the concept applies to three-dimensional as well as to two-dimensional 
flow. The presence of shock waves in a flow of extremely high Mach number can 
change the order of magnitude of the local sound velocity. However, a simple in- 
vestigation shows that this sound velocity cannot be greater in order of magnitude 
than Ma,@, where £ is the inclination of the surface causing the shock and is assumed 
small but of order 1/M or larger. Thus with shocks present, the local Mach number 
will remain of order 6—! or larger. Hence the consideration of the hypersonic flow 
about a slender body as a non-stationary problem in one less dimension remains 
valid when shoek waves and the resultant entropy changes are present. 

The general similitude may be expressed thus: If a slender body of the shape 


g(é, 7, 5) = 0 (8) 
where 
x = bé, (9a) y = dn, (9b) z= 6 (9c) 
is placed in a uniform stream of large Mach number the problem is identical with a 
non-stationary problem in y, z, and ¢ where 
t = x/Mao (10) 


and is characterized by the parameter K as given by Eq. (3). The boundary condi- 


tion satisfied on the surface is 


bs v w 
Kg: Ti —)}% +t(—)s.= 0 (11) 
ao ao 


where v and w are the velocity components in the y and 2g directions, respectively. 


The drag and lift coefficients based upon an area of magnitude bé are given by 


Eqs. (5a) and (5b). 
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ON THE SUPERPOSITION OF A HEAT SOURCE AND 
CONTACT RESISTANCE* 
By S. A. SCHAAF (New York University) 


1. Introduction. The usual condition for a thermal contact resistance! is that the 
temperature discontinuity 7,— 72 across the interface of two heat conductors [1 | 
and [2] should be proportional to the rate of heat transfer H per unit area there, i.e. 


T, = T2 = RH, 


a 


= “resistance” constant. (1) 


If heat is generated at the same interface, however, the quantity H in (1) cannot be 
interpreted unless a more specific physical model of the interface is considered. For a 








{1] \ ¢ [2] 
, \ Uy 
\ 7: 
\ 








Fic. 1. Magnified view of “lubricated friction” model for interface, consisting in a source S between 
the two media and separated from them by contact resistances R; and Re. A typical temperature dis- 
tribution for some fixed time t>0 is superimposed. 


large class of problems, including those in which the heat source is caused by friction, 
appropriate physical models of the interface lead to a condition of the form 


T, —T2= R(CiS = HM) = R(C2S - H2) C,+C, = 1, (2) 


where S is the rate of heat generated per unit area at the interface, H; is the rate of 
heat per unit area flowing into the corresponding medium and C; is a constant whose 
interpretation depends on the nature of the interface. Condition (2) is obtained by 
considering models of the interface in which either the contact resistance or the heat 
source is broken up into two parts. If the source is caused by friction as the two solids 
slide against each other with the interface as slip plane, these two models correspond 
respectively to lubricated and to dry friction. 


* Received May 27, 1946. 
1W. A. Mersman, Heat conduction in an infinite composite solid with an interface resistance, Trans. 
Amer. Math. Soc. 53, 14-24 (1943). 
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Lubricated friction. Here we suppose (see Fig. 1) that the heat source is located 
between two known contact resistances R; and Re, with R,;+R2.=R. This is the case 
if the two solids are separated by a layer of lubricating fluid, for example melted 
material from one of the solids. Heat is then generated in the turbulent fluid and flows 
to the two solids through film contact resistances of the ordinary kind. If we suppose 
that the fluid temperature is 7, then according to (1) we have 


T-—T;= RH; at x=0, (ti = 1, 2). 
Eliminating 7, using Ri +R.,=R and H,+H2=S, we obtain 
T, ~ T2 = RS ad RHA, soe (RS “a RH») at += 0, 


which is equivalent to (2) if we put C;=1—R,/R. 
Dry friction. An alternative model for the interface (see Fig. 2) is to suppose that 


[1] (2) 














os] 
Source Sy 
Resistance R 











Fic. 2. Magnified view of “dry friction” model for interface, consisting in two sources S; and Sz 
separated by a contact resistance R, with typical temperature distribution as some fixed time ¢>0. 


the heat is generated on each of the two surfaces, so that there are two sources 5S;(¢) 
and S(t), with S;+.S,=S, separated by a single contact resistance R. This seems a 
reasonable model for dry friction, where the heat is generated on the two sliding sur- 
faces which are separated at most points by a very thin air gap which constitutes the 
contact resistance. It further seems reasonable to suppose that these two sources are 
proportional, so that 


Si(t) = C;S(2), So(t) = C2S(t), Cy + C2 = # 


where C; is an empirical constant depending on the elastic and cohesive properties of 
the two media and on the geometry of the surface roughness, but not on the source 
strength S(#) nor the temperature. The rate of heat per unit area flowing across the 
resistance from [1] to [2] is then S(t) —#,(t), so that according to (1) we have 


T; —— T> = R(S; wae Hi) = R(H2 aa Se) at x=0@ since Sy + Se = Hy; + Ho. 


This is equivalent to (2) providing that the assumption S;=C;,S is valid. 
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In Sec. 2 we solve a typical heat transfer problem whose formulation involves the 
boundary condition (2). 

2. Problem and solution. Let us consider the case of two semi-infinite media [1] 
and [2], initially at zero temperature and in contact along the plane x =0, where there 
is a heat source S(t) and a contact resistance R. Then, in the usual notation, we have 


the following boundary value problem: 


0°T; 1 OT; , , 
—_ = — —» (¢= 1,2), x>0, t>0. (3) 
Ox? a; dt 
T;=0, 0 ‘2£>@ (4) 
aT, OTs 
ty ~ Toa R[ csi + Ki “| on ~ BCS + Ke— b ; 
Ox Ox (5) 


x= 0, i> 0, C,+C.= 1 


Equation (5) is the same as (2) since 


oT; 
H; = — K; at x=0, t>0. (6) 





Ox 


The solution of (3) and (4) in terms of the undetermined functions H;,(¢) is well 


known? to be 
t . ov 
7; = a: f H(t — u) exp (— yi/4u)du/V u, (7) 
0 7 


where y;=x; ‘s/ ai, Bs=*/ a1) Kw/r. The functions H;(t) may then be determined by 
use of the Laplace transform from the integral equation (5) and the additional re- 
quirement 

H(t) + H2(t) = Sd), (8) 
which is obtained by equating the two right members in (5) and then using (6). We 


denote by f*(s) the transform of f(t), i.e. 


io) 


fr(s) = f e-**f(i)at (9) 


and recall that 
t * 
| f f(t - w)gtwau = f*(s)-g*(s). (10) 
0 


Transforming (7), we obtain 
T# = B:H#*(s)-Vx/s exp (— vir/s), (7)* 
since® 


— 


exp (— y,/4t) i —_ - 
[=e ee) = \/n/s exp (— yivV/S). (11) 
V 


2H. S. Carslaw, Conduction of heat in solids, p. 153. 
3 J. Cossar and A. Erdelyi, Dictionary of Laplace transforms, ACS-DSRE No. 71, p. vi-20. 
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The transforms of equations (5) and (8) are 
Bil #*\/2/s — BoH#\/x/s = R(CiS* — H#), (5)* 
Hy} + H# = S*. (8)* 
Solving for Hi#*, we obtain 
re bi tas SO (12) 
(81 + BV a + RVs 
We will treat the cases S=constant So, and S=S(t) separately. In the former case, 
we have 
S* = So/s, (13) 
and the inverse of (12) is then known‘ to be 
BoSo So(B1C1 — B2C2) 


+ — -exp [(81 + B2)2rt/R?] erfc [(81 + Be) rt/R], (14) 


ae = 
Bi + Bo 


- 81 - Bo 


where 


9 20 
erfc (x) = —f e~“ du. 
VJ: 


H2(t) may be obtained from (14) by interchanging subscripts 1 and 2; the expressions 
for 7, then follow immediately from (7). We observe that as t>«, H; approaches 
the constant value it would have for R=0, i.e. for no contact resistance. It follows 
that, although 

T (0, t) ~ const. y t as to, 
the temperature discontinuity 


T,(0, t) — T2(0, 4) > const. as to. 


For the case S=S(t), we substitute (12) into (7)* to obtain 


Boe + CiRV/s e exp (— yiV/s) 
T* = B,S* EF Bova + OiRV/s Vm exp - yaN | (15) 


(B: + B2)V/a + Rvs Vs 
This is in the form of a product of two functions with known inverses,® so that (10) 


may be used to obtain 


es : C; exp (— yy 4u) m(B2C2 — BiCs) 
T,= a f S(t — “u)- <——————= + - ae 
0 \/ U R 


[Bi + Be = (81 + Be)? faite — v1 ' 
-exp | ——- wiv + —— ou | -erfc ——— ru + baw. (16) 
R , R 2V/u 


The expression for JT; may be obtained from (16) by interchanging the subscripts 1 


4 J. Cossar and A. Erdelyi, loc. cit., p. vi-76. 
5 J. Cossar and A. Erdelyi, loc. cit., pp. vi-77, 78. 
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t/y 


[steel] [copper] 














Fic. 3. Temperatures distribution in steel and copper media, separated by a contact resistance and a 
heat source whose strength is proportional to the time, S=7t. 


and 2. The evaluation of (16) is given graphically in Fig. 3 for the typical numerical 
case (in c.g.s.) B:=1.8 (steel), Be =.6 (copper), Cx =C,2=1/2, R=4.2, S(t)~t; which 
corresponds approximately to the early stages of heat transfer between a steel gun 
barrel and a copper shell under constant acceleration, 
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W. PRAGER 





























